UH-511-1066-05 



Fully differential Higgs boson production and the di-photon signal 
through next-to-next-to-leading order 

Charalampos Anastasiou 

Institute for Theoretical Physics, 
ETH, 8093 Zurich, Switzerland 

Kirill Melnikov 

Department of Physics and Astronomy, University of Hawaii, 
2505 Correa Rd., Honolulu, Hawaii 96822 

Frank Petriello 

Department of Physics, Johns Hopkins University, 
3400 North Charles St., Baltimore, MD 21218 

Abstract 

We describe a calculation of the fully differential cross section for Higgs boson production in 
the gluon fusion channel through next-to-next-to-leading order (NNLO) in perturbative QCD. The 
decay of the Higgs boson into two photons is included. Technical aspects of the computation are 
discussed in detail. The implementation of the calculation into a numerical code, called FEHiP, 
is described. The NNLO i^-factors for completely realistic photon acceptances and isolation cuts, 
including those employed by the ATLAS and CMS collaborations, are computed. We study various 
distributions of the photons from Higgs decay through NNLO. 
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I. INTRODUCTION 



Perturbative calculations in quantum field theory have been performed since the birth 
of QED in the 1930s. Enormous progress has been made since this early era, and the tools 
used to obtain these results have improved dramatically. Much of this progress has been 
spurred by the constantly increasing demands of the high energy experimental program. 
Future experiments such as the LHC and the Next Linear Collider are now demanding 
perturbative calculations to next-to-next-to-leading order (NNLO) in the relevant coupling 
constants. These results must also be directly applicable to experimental measurements, 
which requires calculational algorithms flexible enough to allow arbitrary constraints on the 
final state of reactions. 

For a long time, virtual corrections hindered the progress of perturbative calculations. 
With the advent of computer algebra and the realization that loop integrals satisfy simple 
identities that follow from their generalized hypergeometric nature virtual corrections 
for processes with a small number of kinematic invariants became relatively straightforward 
[1; IE 1^ 1^- However, physical results require the inclusion of real emission corrections, in 
which additional massless partons are emitted into the final state. It was recognized early in 
the history of perturbative calculations, first in QED and then in QCD, that such emissions 
lead to infrared and collinear singularities when the initial or final state configurations 
become degenerate. The Bloch-Nordsieck and Kinoshita-Lee-Nauenberg theorems 0, 0, S| 
show that both virtual corrections and real emissions must be included to obtain a physically 
meaningful result, because divergences only cancel when these components are combined. 
The methods cited above rely upon special features of loop integrals. They can be used to 
obtain virtual corrections, or inclusive cross sections related to virtual corrections through 
the optical theorem, such as the total hadropro duct ion cross section in e~^e~ collisions. 
Recently, these methods have also been extended to phase- space integrals for total cross 
sections and simple kinematic distributions [1, U, 13, IS 3 However, observables 
where the phase space is non-trivially constrained cannot be obtained using these techniques. 
Unfortunately, these are exactly the quantities required by experiment; because of final-state 
cuts, inclusive results are of limited use. 

Virtual corrections possess a simple mathematical structure, and the extraction of their 
singularities proceeds in an observable-independent fashion. With real corrections, different 
kinematic cuts are imposed on the final state for different observables. Since we would like to 
perform calculations valid for arbitrary cuts, the extraction of singularities from real emission 
corrections should be performed in the presence of an unspecified "measurement function" . 
The factorization of soft and collinear emissions renders such an extraction possible. The 
resulting cancellation of singularities requires that the measurement function allows only 
"infrared safe" observables that can be computed in perturbation theory [l6| . 

An efficient algorithm for extracting singularities at NLO in perturbative QCD was con- 
structed in [13, [13, advancing earlier work on the subject [l^. This "dipole-subtraction" 
algorithm identifies universal soft and collinear counterterms, called dipoles, that can be 
subtracted from the real emission contribution to an arbitrary process to make it finite. 
The dipoles can then be analytically integrated over the restricted phase space, since the 
measurement function for infrared safe observables collapses to a measurement function of 
lower multiplicity in the soft and collinear limits. After integration, the dipoles cancel the 
infrared divergences in the virtual corrections, leading to a finite result. 

For processes where the perturbative corrections are large, or control of the theoretical 
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error is crucial, perturbative calculations must be extended to NNLO. Either an extension 
of the dipole formalism to NNLO or the development of an alternative approach is needed 
to enable these computations. Significant effort has been devoted to a generalization of the 
dipole subtraction method ji^ 21, 2^, IS 24, 2^ 2^ 27. lisl. This extension has remained 



elusive. Although the infrared behavior of loop amplitudes |29l . |30| and the infrared limits of 
real emission corrections are universal at NNLO [3l|, it becomes much more difficult to dis- 
entangle the singularities of two unresolved emissions, and to construct process-independent 
counterterms. 

We have recently developed an alternative approach to the problem of real radiation 
at NNLO js^. Our method differs conceptually from the dipole subtraction approach in 
two important ways: the finding of singular phase-space regions is completely automated; 
the cancellation of the 1/e poles which describe divergences in dimensional regularization is 
performed numerically, and no analytic integrations are required. These features guarantee 
that it can be used to extract and cancel singularities at any order in perturbation theory, at 
least in principle. The primary complication that occurs at higher orders is the presence of 
overlapping singularities. These can be disentangled using sector decomposition 3^ S^, 3^ ■ 
Existing symbolic manipulation programs, such as MAPLE or MATHEMATICA, provide 
a suitable framework in which to program the algorithm. The extraction and cancella- 
tion of singularities is thus achieved in a completely automated fashion, with little human 
intervention. 

We have used this approach to perform two fully differential calculations through NNLO 
in QCD: the 0{al) correction to e~^e~ 2 jets 32, 3^, and the Higgs boson hadroproduc- 
tion cross section through gluon fusion ji^. Both calculations permit arbitrary cuts on the 
final states and can easily be extended to include decays of the final state particles. These 
successful applications demonstrate the vitality of our approach and its potential relevance 
for other problems of phenomenological importance; we therefore believe it is important to 
describe it in a simple, pedagogical fashion. We attempt to do so in this manuscript. Our 
discussion will be centered around the calculation of the Higgs hadroproduction cross section 
completed recently 0]. We will extend this result to include the decay of the Higgs through 
the channel H — > 77. We therefore pursue two goals in this paper: a thorough discussion 
of the analytic aspects of our calculation and its implementation into a numerical code, and 
a presentation of phenomenological results for pp^H-\-X—>-'~f'j + X through NNLO in 
QCD. 

This paper is organized as follows. In the next Section we describe phenomenological 
issues relevant for Higgs boson hadroproduction. In Section IIIII we introduce our notation, 
and describe the basic setup of our calculation. In Section |TVl we discuss Higgs production 
in association with up to one parton. Section |3 is devoted to the calculation of the coUinear 
counterterms. The phase-space parameterization is a crucial element of our approach; we 
describe details regarding the choice of parameterizations for Higgs hadroproduction in Sec- 
tion |^I We discuss how to handle the various forms of singularities which appear in NNLO 
computations in the following Section. In Section IVIIII we discuss how the matrix ele- 
ments that appear in the double real emission contribution for Higgs hadroproduction are 
treated in our calculation. We then present phenomenological results for pp ^ H + X and 
pp—i-H + X-^'j'j + X through NNLO in perturbative QCD. We next describe the imple- 
mentation of our results into the numerical code FEHiP. We conclude with a discussion of 
the advantages and disadvantages of our approach, and identify directions for future work. 
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II. HIGGS BOSON PRODUCTION AT HADRON COLLIDERS 



We describe here the aspects of Higgs phenomenology at the LHC needed for our calcu- 
lation. There are several mechanisms for Higgs production at hadron colliders (see jsll for 
a review). For most Higgs masses, gluon fusion through a top quark loop, gg —>■ H, is the 
dominant production mechanism. For Higgs masses in the range preferred by the global fit 
to the precision electroweak data Isll, rrih ~ 110 — 130 GeV, the gluon fusion cross section 
is approximately 60 pb. The Higgs branching fraction into two photons is larger than 10~^ 



in this mass range |40|, |41[. The search strategy |42| for the Higgs signal is then to look for 
events with two isolated photons and reconstruct the mass of the Higgs boson by studying 
their invariant mass distribution. The photons are required to have transverse momenta 
p^l^ > 40 GeV and p^f^ > 25 GeV; they must also be produced in the central rapidity region 
\ri\ < 2.5 lij. The major irreducible background to two photon events is direct (prompt) 
di-photon production in hadronic collisions a 1% photon energy resolution is needed 
to distinguish the if — 77 signal over the background. Typically, an isolation cut is also 
imposed upon the photons; this suppresses photons from the decays of large p± hadrons, 
such as ir^ —>■ 77, and from the poorly known fragmentation production of prompt photons. 
Several possible isolation cuts have been proposed 4^, 4^, 4^ 4^. The simplest possibility 
is to require that a photon candidate does not have additional transverse energy -ET.min de- 



posited within the region i?is = y {v ~ V-yY + (0 ~ ^i)^ the (r^, 0) plane. Typical values 
used in previous studies are -ET.min = 4 — 15 GeV and Ris = 0.4. 

The Higgs production cross section receives large perturbative corrections and depends 
strongly on the renormalization and factorization scales. For example, for rrih = 100 — 
130 GeV, the cross section f or pp ^ H + X at ^/s = 14 TeV increases by a factor 1.5 — 1.7 
when the NLO QCD corrections are included ji^,!!^. The residual scale dependence at NLO 
is approximately thirty percent. This peculiar behavior of the perturbative series motivated 
several NNLO calculations of the inclusive Higgs production cross section jOi. i49i..5(li5Li52i|. 
These studies found no breakdown of the perturbative expansion; while the NNLO effects 
are sizable, they are smaller than the NLO ones. The NNLO cross section is also much more 
stable against variations of the renormalization and the factorization scales. These results 
were confirmed ^] in the framework of threshold resummation, which exploits the fact that 
because of the large value of the gluon density at small Bjorken x, the Higgs production cross 
section at the LHC is dominated by the partonic threshold, i.e., the z —>■ 1 region, where 
z = m^/spart- The terms that are singular in this limit can be systematically resummed 
to all orders in as, and the results compared to the complete NNLO calculation. The two 
approaches agree well, indicating that the uncalculated higher order effects are likely to be 
within the uncertainty assigned to the NNLO result. 

Much is also known about less inclusive quantities for Higgs boson production. The NLO 
p± and rapidity distributions for Higgs boson production at high p± are computed in |5j| . In 



addition, the rapidity distribution of the Higgs boson has been computed through NLO [11 
4^. The p± distribution of the Higgs bjison has been investigated using various resummation 



formalisms by different groups |55|, |5a, |57|, |58[. Monte Carlo event generators accurate 
through NLO for the Higgs + jet process have been published [H^ [H^. Higgs production 
is also included in existing shower Monte Carlo event generators, such as PYTHIA and 
HERWIG, and in the MC@NLO event generator that correctly combines single hard gluon 
emissions with the HERWIG parton shower . The di-photon invariant mass distribution 
in pp collisions, including both the signal pp^H + X—>-'y'~f + X and the prompt photon 
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background pp — > 77, can be found in 0, 4^ Q^. Ref. [i^ presents the partonic level event 
generator DIPHOX, where both direct and fragmentation components of the prompt photon 
production are computed through NLO in perturbative QCD. In |46i], the QCD corrections 
to the gg 77 channel of diphoton production are computed; combined with DIPHOX, this 
presents the most up-to-date analysis of the two photon background to Higgs production 
and enables a careful analysis of the signal-to-background ratio as a function of the isolation 
cuts. In jo^], the interference between the signal and background is shown to be negligible. 

Clearly, a substantial amount is known about Higgs hadroproduction; unfortunately, only 
the inclusive cross section is known through NNLO. Since the NNLO corrections are large, 
it is desirable to also know differential quantities to this order. This is necessary to compute 
the i^'-factor, -ft^NNLO = ctnnlo/clo, that corresponds to realistic experimental acceptances 
with cuts on the photons and jets in the final state. Although the kinematics of the if — >■ 77 
decay is not altered by higher order QCD effects^, the QCD corrections to the production 
process change the kinematics of the produced Higgs boson and lead to modifications in the 
kinematics of the produced photons. In order to compute the relevant acceptance through 
NNLO, the Higgs boson production cross section must be known at the differential level. 
We begin our description of this calculation in the following Section. 



III. NOTATION AND SETUP 

We study the production of a Higgs boson with momentum ph in the collision of two 
hadrons, hi, /i2, carrying momenta Pi,P2'- 



hi{P{) + h2{P2)^H{ph)+X. 



(1) 



Within the framework of QCD factorization, the cross section for this process can be written 
as an integral over hard scattering cross sections cTij for the production of the Higgs boson 
from the quarks and gluons, multiplied by parton densities describing the distribution of 
these partons inside the colliding hadrons: 



/■I 

^ = Ylj dXidx2fl:'^^\xi)ff''\x2)(yij^H+x{Xi,X2). 



(2) 



The sum is over the parton flavors i,j in the hadrons hi,h2, and fi^^\ f'f'^'' are the cor- 



responding parton densities. The initial-state partons z,j for the hard scattering partonic 
process carry momenta pi = XiPi and p2 = X2P2- The partonic cross sections for the pro- 
cesses i + i H + X, can be calculated perturbatively; here, we compute them through 
in the strong coupling expansion: 



a. 



ij^H+X 



(2) 



(0) , Ois (1) , / 



(3) 



The partonic cross sections aij contain divergences arising from intial-state collinear ra- 
diation; these are removed by recasting the parton-densities in the MS-factorization scheme. 



^ Note that this statement is violated by the decay of the Higgs boson into two gluons and two photons, 
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The hadronic cross section of Eq. |21 is computed in terms of finite parton densities / and 
finite partonic cross sections a, 



^ = Y1 j dXidx2fi (Xi)/] -'[X2)Crii^H+X[Xl,X2) 

The finite and "bare" parton densities are related via 



(h) 



(4) 



(5) 



where we have introduced the convolution integral 



{f(^g){x)= dydzf{y)g{z)6{x~yz). 
Jo 



(6) 



The functions Fj, are given in the MS scheme by 



Tij{x) = Sij6{l- x) 



a. 



TT e 




(7) 



where the Altarelli-Parisi kernels P/"^ can be found in |63l. IH^]. We note that the complete 



NNLO corrections to these kernels have recently been computed j65|- e = (4 — (i)/2 is 
the usual dimensional regularization parameter; all calculations in this paper are performed 
using this regularization scheme. Substituting Eq. (jUI) into Eq. (@)) and comparing with 
Eq. (j2D we find 

'^ii = E / dyidy2Tik{yi)Tji{y2)aki{xiyi,X2y2)- 



(8) 



We compute the finite partonic cross sections o"j,- by expanding 



(9) 



and solving Eq. (jH)) in terms of the coefficients a^^^ order-by-order in the strong coupling 
expansion. In this procedure, we need to consider the convolution integrals of the partonic 
cross sections with the Altarelli-Parisi kernels at each order in the perturbative expansion; 
this will be discussed in detail in a later Section. 

We now discuss the Lagrangian which describes Higgs boson production. As mentioned 
in the previous Section, we consider the gluon fusion mechanism for Higgs production. 
The Higgs coupling to two gluons is induced by a top quark loop if there are other 
heavy quark doublets that acquire mass from the Higgs mechanism, they might also give a 
substantial contribution to the effective Hgg coupling. We focus here on a light Standard 
Model Higgs boson whose mass is smaller than twice the mass of the top quark: rrih < 
2mt ~ 350 GeV. The interaction of the Higgs boson with two gluons can then be described 
by a point-like vertex |4^; this is formalized by introducing the effective Lagrangian 



(10) 
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where is the gluon field strength tensor, H is the Higgs field, and v ~ 246 GeV is the 
Higgs boson vacuum expectation value. The Wilson coefficient Ci and the renormalization 
factor Zi, defined in the MS scheme, are 




+ 



TT 











_e2 


e 



2 r2777 19 

288 ^ 16 * 



O a 



67 1 
96^3 * 



(11) 



where = dsil^) is the QCD MS coupling constant, defined in the theory with Uf 
fiavors, and 



log 



m 



A 







top , 



11 



6 ' 



Ai 



51 19 



n 



8 24 



(12) 



The Feynman rules for the ggJi vertex follow from the effective Lagrangian in Eq. |TO|) . The 
Higgs boson couplings to light fermions are neglected. The leading order partonic process 
is then gg H] at higher orders in perturbation theory, other partonic processes also 
contribute. 

At the LHC, the energy of hadronic collisions is much larger than the Higgs mass; it is 
therefore not obvious that the approximation of a point-like Higgs coupling to two gluons 
is sufficiently accurate. However, Higgs bosons with masses mt ~ 100 GeV are predomi- 
nantly produced close to the threshold of the partonic collision, with an average transverse 
momentum of tens of GeV. The kinematic invariants in the partonic gg H + X pro- 
cess never become large enough to resolve the top quark loop, unless the large p± region is 
specifically probed. Since the contribution from this region is negligibly small, the point-like 
approximation is valid. 

We now discuss what partonic contributions are required to compute the differential cross 
section for Higgs hadroproduction at NNLO. At the partonic level, the leading order process 
is gg — > H. At NLO, three other partonic processes appear: gg — >■ Hg, qg Hq and 
qq Hg. At NNLO, we must compute: i) two-loop virtual corrections to gg H; ii) one- 
loop virtual corrections to gg — > Hg, qg — * Hq and qq — > Hg; iii) inelastic processes with 
two partons in the final state: gg — » Hgg, gg Hqq, qg Hqg, qq — > Hgg, qq — > Hqq, 
and qiqj — > Hqiqj. Each of these contributions has the generic form 



a = J dUi |A^ij^/f+Zpartons| Fj{xi, X2,Pl,P2, {?«}), 



(13) 



where dUi denotes the integration over the Higgs phase-space and the phase-space of I 
additional partons in the final space, Aiij includes both the matrix elements and any required 
loop integrations, and Fj describes the observable under consideration. 

The loop integrations are universal for all observables, and can be performed with well 
established methods. We use integration-by-parts identities and recurrence relations 0] to 
calculate the virtual corrections to both the LO and NLO processes. The recurrence relations 
are solved using the algorithm described in P] and implemented in [H^. The resulting 
master integrals are then evaluated directly. This is straightforward for the two-loop virtual 
corrections, since the / = phase-space integration just gives 5(m| — {pi + P2)^)', for the 
virtual corrections to the NLO processes, the resulting master integrals have to be integrated 
over the / = 1 particle phase space. This can produce additional singularities. Care must 
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be taken to assure that all singularities are extracted properly. However, because the two- 
particle phase-space is simple, the extraction of singularities is easy and proceeds along 
the lines discussed in Hence, dealing with either two- or one-loop virtual corrections 
to Higgs hadroproduction is straightforward; we will discuss these briefly in the following 
Section. 

The situation is drastically different for the double real emission channels. As discussed 
in the Introduction, efficient extraction of infrared and collinear singularities from this com- 
ponent is still an open issue. It is known how to compute analytically the phase-integrals for 
the total cross section [1, |5l|, S^], where we must set Fj = 1. It is also possible to compute 



analytically simple kinematic distributions where the measurement function takes a simple 
form. For example, to compute the rapidity Y distribution of the Higgs boson in the frame 
of the two hadrons we must insert 




pI + pI 

Ph - Ph 




X1P2 ■ Ph 
X2P1 ■ Ph 



(14) 



where the z-axis is the beam axis. Phase-space integrations of this type, where the mea- 
surement function can be written as a delta-function constraining a covariant quantity, can 
be mapped to loop integrals and solved using the techniques discussed in the previous para- 
graph [12, 12,]. The rapidity distribution has been computed analytically for electroweak 
gauge boson production at hadron colliders 0, these computations require similar 



phase-space integrations as for Higgs boson production. However, the measurement func- 
tion Fj can take very complicated forms, which are unsuitable for an analytic evaluation of 
the cross section, if additional components of the Higgs boson momentum or the final-state 
partonic momenta are probed, a jet finding algorithm is applied, or the decay of the Higgs 
boson with all relevant experimental cuts is included. 

The difficulties related to the evaluation of the double real emission components can be 
summarized as follows. Naively, Eq. (jl3p is finite for these contributions, and the limit e ^ 
can be taken. However, this is only true for non-exceptional momentum configurations. If 
the momenta of some particles become soft, qi 0, or collinear, ■ qj — > 0, qij 7^ 0, 
the matrix element diverges and can not be integrated in four dimensions. Computing the 
contribution of the real emission graphs to the total cross section amounts to integrating 
Eq.lfT^ over the entire phase space {Fj — >• 1), so that the soft and collinear regions are 
included. However, for the differential cross section such an integration is not allowed, since 
we want to keep the kinematics of the final state intact. The challenge is then to extract 
the singularities from Eq. (jl3p in a way that correctly accommodates both singular and non- 
singular limits, and does not require any integrations to be performed. This can be done 
using the approach suggested in js^], which we explain in detail in this paper. We sketch 
here its outcome. Using this method, we are able to rewrite dureai in the following form: 

a,„,.^M^. (15) 

where Ai are functions non-singular everywhere in the phase space. No specific information 
about the measurement function Fj is used in this derivation. The functions Ai contain no 
residual e-dependence, and can therefore be computed numerically in four dimensions. The 
expression for the double real emission component of Eq. (|15p is then combined with similar 
expressions for the virtual corrections, and the singularities in e are canceled numerically. 
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FIG. 1: Examples of diagrams that contribute to the g + g ^ H cross section 

After the cancellation of singularities is established, we can drop the singular terms from the 
expression for the cross section and implement the finite part into a numerical code. This 
is the basic strategy which we discuss in detail in the remainder of this paper. 

IV. PRODUCTION OF THE HIGGS BOSON IN ASSOCIATION WITH UP TO 
ONE PARTON 



We start with the partonic cross sections for producing the Higgs boson and no partons 
in the final state: 

g{pi)+g{P2)^H{pf,). (16) 
The 2 — > 1 phase-space is simple, because of momentum conservation. We derive 



y dHo = j d'^Ph5\ph - Pi -P2)5{pI - ml) = 6{ml - s). 



(17) 



where s = {pi + P2)'^ is the partonic center of mass energy squared. 

The most complicated part in computing the partonic channel Ugg^u is the evaluation 
of the virtual corrections through two- loops (see Fig. Q). Fortunately, these corrections are 
known from the arialytic calculation of the inclusive Higgs boson production cross section 
through NNLO 0, S H , and we use these results in this paper. 

We next study the cross sections for partonic processes with the Higgs boson and a quark 
or a gluon in the final state: gg Hg, qg — > Hq, and qq — >• Hg. For these processes, 
we must compute the corresponding tree-level and one-loop amplitudes. We consider the 
process g{pi) + g{p2) H{ph) + g{pz) as an example. Typical diagrams are shown in Fig|21 

Consider a contribution arising from the interference of two tree-level diagrams to the 
differential cross section. It can be written as: 



Ar(si3, S23, i^j) 



'S13S23 



(18) 



where Sij = {pi — pj^ and Fj is the measurement function which defines the observable we 
want to compute. The only information we need about the numerator in Eq. ()18|l is that 
it is a finite function in the limits S13 — > and S23 — 0. The structure of the infrared 
and collinear singularities is fully determined by the denominator of Eq. ()18|l . We use this 
observation to derive an expansion in e for this denominator in terms of delta functions 
and plus distributions. Having done that, we treat arbitrary numerators using a numerical 
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FIG. 2: Examples of diagrams that contribute to the production of the Higgs boson in association 
with one parton. 



subroutine for Fj, and defining procedures to compute the action of delta functions and plus 
distributions on integrable functions. 

Therefore, the basic integral we have to consider is 



Igg^Hg = J d'^Phd%5^ 



Ph 



m. 



(5+ (^3) (5^ {pi +P2-Ph- Ps) 



Fj{Si3,S23) 

S13S23 



(19) 



This integral is potentially singular for S13, S23 = 0. To extract the singularities, we pa- 
rameterize the phase-space in terms of variables that range from to 1 in such a way that 
the singularities are mapped to the boundaries of the integration region. A convenient 
parameterization is in terms of the variables Ai, A2, where 



Sis 



ml 1 + Ai 



;i-Ai; 



S23 



ml 1 + A2 



Ai Ai + A2 A2 Ai + A2 

In this parameterization the integral in Eq. ()19|) becomes 



A2). 



(20) 



2s 



J^XidX^d |^AiA2-^^ (1-Ai) 



)-^-^(l-A2) 



-l-e 



1 + A1A2 



(1 + Ai)(l + A2) 



ml 



:i + A0(l + A2) 



AiA2(Ai -f A2) 



^j(Sl3,S23). 



(21) 



The delta function appears because of the momentum conservation S12 + S13 + S23 = ml, and 



prevents Ai^2 from reaching 0. The singularities that occur as Ai, A2 — > 
form. To extract them, we rewrite the singular terms (1 — Ai)^^^'^, (1 



A 



-l+e 



6{X) 



+ E 

n=0 



e 

n! 



log(A)' 
A 



1 are in a factorized 
- A2)^^^'' using 

(22) 



and expand in e. The result contains delta functions and plus distributions and can be 
integrated numerically with the functions Fj. 

Finally, we must discuss the computation of the interference terms of the one-loop and 
tree-level amplitudes for Higgs boson production in association with a single parton. These 
terms require the calculation of the one-loop amplitude, and integrations over the 2 — > 2 
phase-space variables. Using standard reduction methods, we write the one-loop amplitude 
in terms of master integrals that are known analytically. The integration over the phase- 
space then proceeds in a way described above. 
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V. COMPUTATION OF THE COLLINEAR SUBTRACTION TERMS 



In this Section we discuss the computation of the colhnear subtraction terms. We remind 
the reader that, by comparing the hadropro duct ion cross section written through "bare" 
and renormahzed quantities, a relation between the bare (a) and the renormahzed [a) 
partonic cross sections can be derived, as shown in Eq.(jH)). Having computed the bare cross 
section, we then derive the renormahzed partonic cross section by solving Eq.(jH)) order- 
by-order in the expansion in as- We must consider the convolution integrals of the partonic 
cross sections with the Altarelli-Parisi kernels at each order in the perturbative expansion. 
Technical aspects of this computation are discussed in this Section. 

Through NNLO, we have to consider two distinct cases when aij is either the leading 
order or the next-to-leading order cross section discussed in the previous Section. Consider 
first the leading order cross section. Since a'^^'J ~ S{1 — z), where z = m1/s, the required 
convolution integrals are of the form: 



/' dy^ r dy2Ta{yi)T,{y2)6{l - ^) = z (F, ® T,) (z). (23) 

J- ym 



We therefore need to consider convolutions of the splitting functions which appear in the 
perturbative expansion of F^, F;,. The splitting functions generically contain delta-functions, 
plus-distributions and regular functions. The convolutions that involve delta-functions are 
straightforward. The convolutions of two regular functions and the convolutions of a plus- 
distribution with a regular function are performed numerically. The convolution of two 
plus-distributions requires additional analytic work. 
Consider the convolution of two plus-distributions: 



nm 



"ln"(l -x)' 




"ln"^(l-x)' 


1 — X 


+ 


1 — X 



(24) 



We first apply the definition of plus-distributions in the convolution integral: 

Imn = / dydz^ ^— ^ ^ {5{x - yz) - 5{x - y) - 5{x - z) + 5{x - 1)} (25) 

Jo 1 — 1/ I — z 

Eq.((2H) can not be integrated term by term because of divergences; to make such an inte- 
gration possible, we introduce auxiliary regularizations: 

= lim hm — ---- / dydz{l-yr^-il-zr^'^ 



e^o a,b->i d"-a d'H 6"+™ Jo 
{S{x - yz) - 6{x -y)- 6{x - z) + 6{x - 1)} . (26) 



The four integrals in Eq. ()26|) can now be computed separately. The first term 

h = C dydz{l - i/)-^+"^(l - z)-^^'^5{x - zy), (27) 
Jo 

after integrating over y and performing the change of variables z = x + [1 — x)X, takes the 
form 

Ia = {l- x)-^+("+'')^ dX b + (1 - x)A]-"' A-^+'^^(l - A)-^+''^ (28) 
Jo 
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We then rewrite Eq.()28|) using Eq.(j22l). The remaining three terms in Eq. ()26|) are easy to 
compute. After differentiating the result with respect to a,b and taking the hmit a,fe — > 1, 
we expand in e. The leading term in the e expansion yields the convolution of the required 
plus distributions. This solves the problem of computing the collinear factorization terms 
when the leading order cross section 0"^°-' ~ 5(1 — is involved. 

We now proceed to the discussion of the integrals that involve the NLO cross sections 
for the production of the Higgs boson in association with one additional parton. Inserting 
Eq. ipTj) into Eq.® we produce integrals of the form 

dXidX2dyidy2T,iyi)r2{y2)f{XiX2)S (^^/ii/sAiAs - , (29) 

where /(Ai, A2) is the integrand of Eq. ()2H) without the delta function and after the expansion 
in e. The parameterization of the phase-space in terms of the variables Ai, A2 is convenient 
for rewriting the integral Eq. ()29|) through successive conventional convolutions, 

y2A4)ri(y2)/(A3, A4). 



1 / 2\ 1 1 

J d\id\2S ( A1A2 -—] J d\3dyiS{\i - yi\3)Ti{yi) J d\idyi5{\2 
^ ^ 



We note that since we work through the relative order 0{a1), one of the two kernels T{yi^2) 
in Eq. ()29|) is a delta function, so that the corresponding integration can easily be performed. 
After that, we are left with a convolution integral that can be treated along the lines de- 
scribed at the beginning of this Section. It can then be directly used in the numerical 
integration with parton distribution functions. 



VI. PHASE SPACE PARAMETERIZATIONS FOR DOUBLE REAL EMISSION 
PROCESSES 

We can now study the contributions to the NNLO cross section from processes with 
double real emissions. These processes appear only at tree-level. However, the integrations 
over the 2 — 3 phase-space are involved. Our aim in the following sections is to present a 
detailed description of their treatment. 

The first step in using the method described in js^l is to choose a parameterization of 
the double real emission phase space in which the integration region is the unit hypercube. 
This is required in order to use an expansion in plus distributions to extract singularities. 
In principle, any parameterization that accomplishes this is acceptable. In practice, finding 
a convenient one that reduces the number of sector decompositions is important for the 
efficiency of the approach. We will discuss here how to choose a parameterization suitable 
for the topologies which contribute to the Higgs production process. 

For the double real emission corrections to Higgs production at NNLO, we must param- 
eterize a 2 — 3 particle phase space, with one massive final-state particle. We consider 
here g{pi) + g{p2) — ^ H{ph) + gipz) + g{j>i) as a prototypical partonic process, although 
the formulae we derive are valid for all such partonic processes. We consider a fixed energy 
for the partonic collision, (pi + P2Y — ^- The scalar products that appear in the matrix 
elements are 

• Sif = {pi —pfY, where i = 1,2 and / = 3,4; 
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• Sih = {pi - Phf, where i = 1, 2; 

• Shf = {pf + Phf, where / = 3, 4. 

The invariant masses of the form s^f are bounded from below s^f > ml, and therefore do 
not lead to singularities. The 2 — * 3 phase space which we must map to the unit hypercube 
is 

dn = 1 [dps] [dp^] [dph]6'^''^ {pi +P2-Ph-P3- Pa)- (30) 

where [dq\ = dg'^~^/(2go)- We will set the overall scale s = 1 in the discussions that follow; 
it can be restored by dimensional arguments. We also denote a generic invariant mass by 
Sab', the subscripts i, f, h will be used as in the above list. Therefore, 2 = 1,2 and / = 3, 4. 

We denote the variables that describe the unit hypercube by A^; the partonic phase 
space is spanned by four independent variables, so i = 1, ... ,4. The singularity structure 
is dictated by the invariant masses that appear in the denominator of the matrix elements. 
In terms of these variables, the invariant masses can take one of the following three generic 
forms, ordered from most to least desirable: 



a factorized form, in which 



an entangled form, in which 
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(31) 



Sab A1A2 . . . 
1 1 



Sab (Al + A2) 

a "line" singularity, in which 

1 1 



(32) 



IX XI- (33) 

Sab — A2 . . . 



The singularities as Ai,A2 ^ in the factorized form can be immediately extracted using 
the expansion in plus distributions of Eq. ()22|) . We will discuss the remaining two structures 
in detail later in the text, and we only briefly describe here how they are dealt with. In order 
to extract the singularity from the entangled form in Eq. (j32|) . we must sector decompose 
in the variables Ai and A2; this involves splitting the integration region into two sectors, 
Al > A2 and A2 > Ai, remapping the integration limits to [0,1], and then expanding the 
resulting expressions in plus distributions. For the line singularity, in which the singular 
region is along an entire edge of phase space rather than just a point, we must first perform 
an additional variable change to remap the singular region to a point, and then typically 
perform a series of sector decompositions. It is advantageous to express as many of the 
singular structures in a factorized form as possible; this form gives only one sector rather 
than several, which leads to smaller and typically simpler expressions. It also preserves the 
kinematics of the parameterization, i.e., there is only one mapping between the Sij and the 
Aj. In the other forms, each sector has a different kinematics. 

In order to choose a convenient set of parameterizations, we must first study the double 
real emission diagrams that contribute to Higgs production. There are three generic types 
that must be considered: 
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1. Diagrams (Ci) in which the Higgs is emitted from either an internal or final-state line, 
and the singular invariant masses that can appear in the denominator are Sj/; 

1 'ummmm 3 1 ^uG^um^sm 3 

\ 





h 



2. Diagrams (C2) for which the potential singular denominators are S34 and Sj^; 

1 ^Tjmmm^f^ 3 




2 ^ 

3. Diagrams (C3) in which the Higgs is emitted from an initial-state line, and the singular 
invariant masses that can appear in the denominator are Sif and Sih, 

1 mmmm3y 3 




The matrix elements consist of interferences between these three classes of diagrams. We 
can not find a phase-space parameterization in which all of the singular scalar products take 
a factorized form, so it is useful to consider parameterizations tailored to each structure. 
We will discuss two different parameterizations of the partonic phase space, which we refer 
to as the "energy" and "rapidity" parameterizations. 

In the energy parameterization, the invariant masses Sif take on simple, factorized forms; 
it is therefore useful for interferences within the first class of diagrams listed above. It is 
named for the observation that in terms of energies and angles with respect to the initial 
beam axis, Sif = —2EiEf[l — cos^jj]; the soft and coUinear singularities in these scalar 
products are therefore factorized. To derive the expression for the phase space in the energy 
parameterization, we begin with Eq. (j3U|) . and use the momentum-conserving (5-function to 
remove the [dh] integration; we obtain 



dUE = ^J dEsdE^dQsdQ^ E^'^E^-^S [1 



Z + S13 + Si4 + S23 + S24 + S34 



(34) 



where z = m\/ s. It is convenient to continue the calculation in the partonic center-of-mass 
(CM) frame. We introduce the CM frame parameterization 



(35) 



Pi = , P2 = 

P3 = E3 (I,sin6'3,0,cos6'3^ 



2(1.0.-1 

, p4 = i?4 (1, sin6'4cos0, sin6'4sin0, cos^4) 



In the expressions for ps and p^ we have suppressed the additional e-dimensional components 
of the momenta; they can be chosen to vanish in this frame. We use the (5-function to remove 
the E^ integration, which sets 

l-z-2Es 



E, 



2 [1 - Es{l - fii ■ n2)] 



(36) 
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and gives the jacobian 2 [1 — i?3(l — ni ■ ^2)]; in these expressions, 

ni ■ = cos^3Cos^4 + sin6'3sin6'4cos0. 



(37) 



We must now consider the angular integrations dQ^ and dQ^. In terms of the polar angles 
^3 and 64, and the azimuthal angle 0, they become 



dn-, 



dcos6'3 



sin^^. 



dQ^ = dcos6'4dcos0 



n 



d-2, 



sm 



-e-1/2 



d-3, 



where is the solid angle in ci-dimensions: 



T{d/2y 

Using these expressions, we can write the phase space in the following form: 



(38) 



(39) 



dllf 



sin26'4 



dE3dcose3dcos^4cicos</) E^'^Ef-^ 

e-1/2 



sm 



sm 



/ [1 - ^3(1 - n, ■ na)] . 



(40) 



The angular integrations clearly range from [—1, 1]. To derive the limits of the E^ integration, 
we note that E3, E^ > 0; using Eq. ljHBj) . we find that this implies 



< E3 < 



(41) 



We are now ready to map the integration region into the unit hj^ercube. Performing the 
variable changes 



E3 = A1E3+, COS^3 = 
COS^4 = — I + 2A3, COS</ 



-1 +2A2, 



-1 + 2A4, (42) 
we derive the final expression for the partonic phase space in the energy parameterization: 



dUE = nJ dAidA2dA3dA4[Ai(l - Ai)]i-2^[A2(1 - A2)]-^[A3(1 - A3 



(43) 



with 



N 
D 

ni ■ n2 



3-4e /o4+2e 



l-{l-z)\^{l-n^■n2)/2, 

A2 + A3 - 2A2A3 + 2(1 - 2A4)v^A2(l- A2)A3(1-A3) 



(44) 



We note that the factor D which appears in the phase space is bounded from below hj D > z, 
and does not contribute to the singularity structure. All limits Aj — >• 0, 1 are separated and 
regulated by e, as required for the extraction of singularities shown in Eq. 
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We now discuss the singularity structure of the scalar products in this paramterization. 
As desired, the invariant masses Sif take on simple, factorized forms: 

si3 = -(1 - z)Xi{l - A2), S23 = -(1 - z)\i\2, 

si4 = -(1-^)(1-Ai)(l-A3)//^, S24 = -(l-^)(l-Ai)A3/I^. (45) 

The other invariant masses have more complex singular structures. We find 

S34 = (1 - ^)'Ai(l - Ai) (1 - ni • n2) /2/D, 

= -{l-z) {A1A2 + A3(l - Ai) - (1 - ^)Ai [1 - Ai(l - A2)] (1 - • n^) /2} /D, 
S2h = -(1 - z) {Ai(l - A2) + (1 - A3)(l - Ai) - (1 - ^)Ai [1 - A1A2] (1 - ni ■ n^) /2} /D. (46) 

The (1 — ■ factor in the numerator of S34 leads to a line singularity along the edge of 
phase space where A4 = 1 and A2 = A3; setting A4 = 1 in Eq. ljl^ . we can write 



(1 - ni ■ n2) |a4=i 



A2(l-A3)- VA3(1-A2) 



(47) 



The invariant masses sih and S2h do not contain line singularities, but do contain several 
entangled singularities; for example, sih vanishes at the following phase-space boundaries: 
(1) Ai = 1 and A2 = 0; (2) Ai = and A3 = 0; (3) A2 = and A3 = 0. These problems indicate 
that the energy parameterization is not well suited for interferences between diagrams of the 
second and third classes listed above. 

For interferences between diagrams in the second and third classes, it is usually better to 
use the rapidity parameterization. With this choice, the invariant masses S34, Sih, and two 
of the Sif take on simple, factorized forms. This parameterization is the fully differential 
extension of the one used in to calculate the NNLO corrections to Drell-Yan production 
of lepton pairs. 

We begin by making explicit what variables we use to describe the phase space. The 
four independent variables we choose are the invariant masses S34, S13, S23, and the variable 
u = pi ■ Ph/P2 ■ Ph- Anticipating the translation of our partonic expressions into hadronic 
results, we note that u is related to the lab-frame rapidity of the Higgs boson hy u = 
Xie~'^^ / X2i where Y is the rapidity and xi and X2 are the standard Bjorken x variables for 
each initial-state parton. We rewrite the phase space in Eq. ljnUj) as 

dllR = j ds34 dw dsi3 ds23 j[dp'i][dpi][dh]5{s2,i-2p?,- Pi) 5{u - pi- ph/p2- Ph) 

x5(si3 + 2pi-p^)5{s2^ + 2p2-p:i)5^'^\pi+p2-ph-p?.-Pi)- (48) 

It is useful to view the production of the final-state particles ^3, ^4, and ph as an iterative 
process; first, ph and the massive "particle" Q34 = P3 -l-p4 are produced; then, Q34 decays 
into p3 and p^. This motivates the following nested decomposition of the phase space: 

dllR = j (is34 (iw (isi3 (is23 dnidn2, (49) 

with 

dill = J[dh]d'^Q34,6{Qli- S34.)6{u-pi-ph/p2-Ph)S^'^\pi+P2- Ph- Q34.), 

dU2 = J [dps] [dp,]Sisi3 + 2pi ■ ps) S{S23 + 2p2 ■ Ps) ^'^^KQm ' P3 ' Pa)- (50) 
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We evaluate dlli in the CM frame of the pi + p2 system. The ^-functions remove all 
integrations except for those describing the azimuthal angles, which give only an overall 
solid angle factor. We arrive at 



dHi 



Qd-2il + z - s34) 



4fl 



u 



us 



34 



2u{l + z)su + {u- z){l 



uz] 



'l+u) 



(51) 



The e-dependent factor will be needed to regulate singularities in m, and S34, as shown 
in Eq. ()22|) : however, the singular structures of these variables are entangled, and must be 
separated. To do so, note that the bracketed term in Eq.()51|) can be written as 



^(44 - g34)(g34 - ^34) 

ii + uY 



(52) 



with = (r ± t)(l ± rt)/r, r = -y/u, and t = y/z. Eq. lfS^ is the of the Higgs, and must 
be positive definite. We must therefore demand < S34 < S34. We set S34 = Aisj^4, where 
Ai is in the range [0,1], and derive 



dHi 



d-2 



'l + z) 



4(1 + m)2 



1 - 



r(l + z) 



{u -z){l- uz){l - Ai)(l - XiKjKp 

{i + uy 



(53) 



where K. 



p,m 



{r ± t){l ± rt). The factor Km/Kp < 1, and the expression in the curly 
brackets does not vanish, so the Ai — > 1 limit has been separated from u, z. However, the 
limits u — >■ z,l/z and z ^ 1 have not been separated, as is clear from the {u — z) and 
(1 — uz) factors in Eq.()53|): these will appear in the denominator and lead to singularities, 
so they must be dealt with. We first note that to keep p\ > we must have z < u < 1/z. 
We then follow 1^ 6^ and set 

with y in the range [0, 1]. The phase space becomes 



4(1 + m) 



1 - 



r(l + z) 



y{l - y){l - zYil - Ai)(l - XiKjKp) 



(55) 



the limits u z,l/z have been separated from z — >■ 1 and moved to ?/ ^ 0, 1. 

Having discussed dHi, we now consider dH2. It is again convenient to work in the CM 
frame of the pi + p2 system. The momenta can be written as 



p^ = ^(1,0,1), P2 = ^(l,0,-1 
P3 = {E,p±sm9,p±cos9,pz) , 



(56) 



where we have again suppressed the e-dimensional components of the momenta. The S- 
functions present in dH2 constrain E, pz, p±, and 6. Using these to remove the integrations, 
we arrive at 



dH. 



d-3 



1+u) 



uis 



34 



34 



"S34J 



p^sin^O 



(57) 
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To express p±sm6 through the variables S34, S23, S13, and u, we must first define the auxihary 
variables 



Do 



'23 



-934(1 + U)[S23{1 + U) + + U - z] - MS23(g34 - g34)(-S34 " ^84) ± \fTh. 
-16u(l + U)^S23S34(S34 " -534) (^34 - S34)(s23 + S23), 



1/(1 -z)l + 



Ai(l-rt) 



r(r + 1) 

In terms of these variables, we have 



p^sin^6' 



- s 



23 J 



l + u)2(s+3 + Si3)(si3 + Si3) 



4m(44 - ■S34)(S34 - S34) 



(58) 



(59) 



There are several consistency conditions that we must apply to these expressions; these will 
give us the integration limits for S13 and S23. First, to ensure the reality of s^^ in Eq. lj^ . 
we must demand D2 > 0; this implies > S23 > —5^3. Finally, since p5_siii^(6') > 0, it is 
clear from Eg . (1591) that we must require — sj^ > S13 > — 5^3. With these limits, we can map 
the S23 and S13 integrations to the region [0,1]. We set 



■^13/ '^13' 



■S23 — ~'^2S23, ■5i3 — — A4(s^3 

and derive the final expression for dn2 in terms of Ai, A2, A4, and y: 



(60) 



dn. 



^d-32 



-4-2e 



1 - z)Jy{\ -y){\- Ai)(l - KKJK^ 



l/(l-y)(l-z)^AiA2(l-A2)A4(l-A4) 
r(r + t)(l + rt) 



-e-1/2 



(61) 

Finally, we must combine dlli and dn2 as shown in Eq. p9|) . and include the jacobian 
for the transformation / ds^A, du dsi3 ds23 I dXi d\2 dX^ dy. This is simply done using the 
variable changes given above, and we derive the following final expression for the phase space 
in the rapidity parameterization: 



dHR = dAidA2dA3dA4 [(1 - Ai)(l - X^KJK,)^ [XMl - A2)]" 



x[A3(l-A3)]'-'^ [A4(l-A4)] 



-e-1/2 



Kpr/{l + u) 



-l+e 



r(l + z) 



(62) 



We have changed ?/ ^ A3 in this equation for consistency of notation. All the limits Aj ^ 0, 1 
have been separated. The singular scalar products in this parameterization are 

sih = -A3(l - ^) [1 - Air(l - rt)/(r + t)] , 
S2h = -(1 - A3)(l - ^) [1 - Ai(r - t)/r/(l + rt)] , 
S23 = -A2A3(1 - ^) [1 + Ai(l - rt)/r/(r + t)] , 
^^24 = -(1 - A2)A3(1 - ^) [1 + Ai(l - rt)/r/(r + t)] , 
S34 = AiA3(l-A3)(l-^)'(l+M)Virp/r, 
(1-A3)(1-^) 



■513 



Kpr [l + Ai(l-rt)/r/(r + t)] 



Al+A2 + 2(2A4-l)V^1^2 



(63) 
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where 

= Ai(l - A2)(l + uf, A2 = A2(l - \i)r{K, - X,K^). (64) 

We have not included S14 in this hst; it can be derived from the other invariant masses (see 
Eg ■ (j34j) ) . and we will see later that the parameterization can always be chosen so that it 
never appears in the denominator. As claimed above, the invariant masses Sih, S2h, S23, 
S24, and S34 are all in factorized forms; none of the expressions in square brackets for these 
invariant masses in Eq. ()63|) vanish. However, S13 clearly has a line singularity when A4 = 
and Ai = A2; in terms of the Aj, 

so this singularity occurs in the physical region. 

Before concluding our presentation of the phase space parameterizations, and beginning 
our discussion of how to handle the entangled scalar products and line singularities, we 
must discuss what happens to the variable z when we use our results to derive the hadronic 
cross section. When we convolute the partonic cross sections with the parton distribution 
functions to form the hadronic cross section, the variable z scales as ^ — *• m1/{xiX2Shad), 
where Shad is the hadronic center-of-mass energy squared, and the Xj are the fractions of the 
hadronic momenta carried into the hard scattering process. It is clear from the invariant 
masses in Eqs. ()45l46l63p that the matrix elements will contain singularities as ;z ^ 1. How- 
ever, it can also be seen from these equations that this singularity is always in a factorized 
form. The normalization factor N in Eq. ()44|1 contains the factor {1 — z)~'^'^, which regulates 
this limit for both parameterizations. Singularities in z and Aj can therefore be treated 
identically. 



VII. FACTORIZING SINGULARITIES 

After choosing a phase-space parameterization for a given term from the matrix ele- 
ments, we must extract its singularities without actually integrating over the Aj. We then 
have differential distributions in the Aj giving us complete control over the kinematics and 
allowing us to compute arbitrary differential observables. For factorized singularities, this 
is simple; the expansion in plus distributions shown in Eq.(j2H) extracts the singularity as a 
1/e pole multiplied by a (5-function which restricts the integration to the singular region of 
phase smce. The singularities can be cancelled numerically and then discarded, as shown 
in js^, 3^ 3^, leaving a finite, fully differential cross section. Our goal will be to reduce all 



other singularities to a factorized form. 

Before discussing the detailed procedure we use to handle entangled and line singularities, 
we must first explain sector decomposition [s^ 13 113 • This is our primary technique of 
separating and factorizing entangled singularities. We can illustrate the main features of 
sector decomposition with a simple example. Consider the integral 

dxdy- — — ^. (66) 

If we naively apply the expansion of Eq. (j22j) . we would conclude that we can simply Taylor- 
expand the numerator of the integrand, and would arrive at 

/naive = /' dxdy ] {1 + 0{e)} . (67) 
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Integrating the 0{e^) term over x, we obtain 

/naive = /'dy— ^ + 0(e). (68) 

Jo y{l + y) 

This is divergent as ?/ — 0; we have clearly missed a singularity. A simple analysis of Eq. ()66p 
shows that the integral diverges logarithmically in the limit x ^ y —>■ 0; if we attempt to 
study these limits separately, as we did in /naive, we miss this singularity. To deal with the 
singular phase-space region x ~ y ~ we use sector decomposition. 

To introduce this technique, we will consider the same integral as in Eq. (j66j) . but with an 
additional function in the integrand: Fj[sab{x,y)], a measurement function, which describes 
kinematic features of the process such as dependence on the parton distribution functions 
and phase-space constraints. This will allow us to discuss issues that arise when sector 
decomposition is used in realistic calculations. Sector decomposition proceeds in a series of 
simple steps. 

1. Split the integration region into two sub-regions; in the first one, x > y, and in 
the second, y > x. The integral / is written accordingly as the sum of two terms, 
I = Ii + I2. We have 

h= dx dy- — ■ — r^Fj[sab(x,y)], l2= dy dx- — ■ — r^Fj[sab(x,?/)]. (69) 
Jo Jo X + y r "'0 Jo x + y r 



2. Remap each integration to the unit hypercube. In Ji, make the change y' = y/x; in 
/2, set X = x/y. Performing these changes of variables (and rewriting x —>■ x,y —>■ y 
for notational ease), we obtain 

h= dxdy-— — —Fj[sah{,x,xy)], h= dxdy— — —Fj[sab{,xy,y)]. (70) 
Jo [1 + yY Jo (1 + x)^ 

3. The singularities in Ii and I2 are now in a factorized form, and can be extracted with 
an expansion in plus distributions. 

4. If a singularity appears in the x ^ 1 limit, as in 

/= /'dxdyii^^, (71) 
Jo (1 — X + yY 

make the variable change x' = 1 — x to map the singularity to x' = 0, and then apply 
the same steps as above. 

There are several features of this process that should be noted. It is very simple to 
program a computer to perform this routine. There are three operations that are performed: 
first, a search for locations of possible singularities of the integrand; second, a variable 
substitution of the form y xy, as in Ji; third, a factorization such as (x + xy)"^ = 
x'^{l + yY in /i in order to find the overall power of x. These operations can be simply 
performed using symbolic manipulation programs such as MAPLE or MATHEMATICA, as 
can the substitution of Eq. ()22|) needed for the plus distribution expansion. Another feature 
to notice is that the mappings between the variables Ai and the invariant masses Sab are 
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different in each sector, as denoted by the different measurement functions Fj[sab{x,xy)] 
and Fj[sab{xy,y)] in Eq.((7n|. If we build an event generator according to the probabihty 
distribution in Aj, we must account for this different kinematics in each sector. Finally, 
applying this technique increases the expression size each time a decomposition is performed, 
so it is best to limit the number of sectors by choosing a phase-space parameterization that 
factorizes as many singularities as possible. 

There are several other features that we include in our computer routine for sector de- 
composition. 

• It is useful to rotate the external momenta in order to change terms in the matrix 
element to a factorized form. For example, in the partonic channel g{pi) + g{p2) ~^ 
H{ph) + giPs) + giPi), we can perform the rotations pi p2 and ps ^ p^. Con- 
sider a term in the matrix element of the form Fj{pi,p2,P3,p4,Ph)/si3/sih, where 
we have described the measurement function arguments using the external momenta. 
In the rapidity parameterization, this term has a line singularity arising from l/si3, 
as discussed below Eq. ljU^ . However, under the rotation pi <-> p2, it becomes 
Fj{p2,Pi,P3,P4:,Ph)/ S23/ S2h, which is in a factorized form. As long as we account 
for the rotation in Fj, which describes the kinematics, this is permissible. 

• When an integral requires sector decomposition, and is singular in both limits A ^ 0, 1, 
we separate these limits by splitting the integration into the two regions [0, 1/2] and 
[1/2, 1]. We then remap each region to the range [0, 1]. We find that this decreases 
the analytical complexity of the result, and improves the numerical precision. As an 
example, consider the integral 

We will evaluate this integral in two ways: (1) by splitting the x-integration as de- 
scribed above, and (2) by directly applying the algorithm of sector decomposition. 
For simplicity, we suppress the measurement function in this example. Using the first 
method, we split the x-integration and derive / = I^^^ -|- with 

where we have changed x — 1 — x in J*^^-' . The first integral requires a sector decompo- 
sition in the limit x — > 0,?/ ^ 0, while the second one is already in a factorized form; 
we obtain / = 4°) + J^°) + J^^), with 

If = / dxdi/2i-^^^— , /(o)= / dxdy2i-^^^-L|^ . 

Jo ^ (2 + x)2 ' y Jo ^ (l + 2y)2 

(74) 

All singularities have been separated, and the expressions can be expanded in e. The 
important point to notice is that except for the terms that must be expanded in 
distributions (y"^"^^^ in /^°\ x~^+^'^ in and x~^~^^ in J*-^-*), all other terms in the 
integrands are finite throughout the entire (x, y) plane. If we instead directly apply the 
sector decomposition algorithm to / without first splitting the x-integration region, 
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we obtain I = + ly, with 

In addition to the components that must be expanded in distributions, the term 
(1 — x?/)"^^'^ in Ij. is singular in the hmit x,y ^ 1. Although this type of singu- 
larity is integrable, we find that it can lead to numerical instabilities when combined 
with parton distribution functions and phase-space constraints, particularly when this 
region of the integration contributes strongly to the result. It is best to avoid these 
terms using the split described above. We note that when we use this split, we can 
choose to map a; — > 1 singularities to a; — >• using the variable change x = 2(1 — x) 
in the [1/2, 1] region. 

Now that we have discussed in detail all the elements that enter our sector decomposition 
routine, we can present our algorithm for factorizing entangled singularities for a given term. 

1. First, check to see if the term can be rotated into a factorized form. This can be done 
by establishing a priority for each denominator structure. In the example given above, 
1 / S22,l S2h would be given a high priority, while l/sis/si/^ would be given a low priority. 
Perform all possible rotations on a given term, and check to see if any of them give 
a high priority integral. If so, then we have a factorized form, and we can expand in 
distributions as in Eq. ()22|) . and go to the next term. 

2. Split each A, integration into the two ranges [0, 1/2] and [1/2, 1], in order to separate 
Aj — > 0, 1 singularities. This can either be done for all Aj automatically, or for a given 
denominator structure we can only split integration regions for those Aj that have both 
Aj and Aj — 1 singular limits. We will assume in the remainder of the discussion 
that the first has been done, so that all singularities occur as Aj — > 0. 

3. Pick two A variables, e.g., Ai and A2, and check whether the term needs sector decom- 
position in these two variables. A term needs sector decomposition if the following 
conditions are met: (1) it doesn't vanish as Ai — *■ 0; (2) it doesn't vanish as A2 — > 0; 
(3) it vanishes as both Ai,A2 0. This is a simple check to program in symbolic 
manipulation packages. It follows from the expressions for the invariant masses given 
in Section II that we can restrict ourselves to singular regions where two of the Aj 
variables vanish; there is no need to consider triple or quadruple singular limits. 

4. If a term doesn't need sector decomposition, return to step 3 and pick another pair of 
Aj. If it does, then apply the technique discussed above. This gives two sectors. Pick 
the first sector, and begin at step 3 for this term. 

5. Eventually, a term will have no pair of variables that requires sector decomposition. 
This term is then in a factorized form, and we can extract phase-space singularities 
by expanding in plus distributions. 

6. Repeat steps 1-5 for all terms. 

We must now discuss how to handle line singularities. These arise when a singularity 
is mapped to an edge of phase space rather than just a point because of a specific choice 
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of parameterization. Our basic method will be to reduce them to an entangled form, and 
then apply the method detailed above. These singularities appear in the more complicated 
topologies, with a larger number of invariant masses in the denominator. In simpler topolo- 
gies, such as the example discussed above, they can be rotated away. It is difficult 
to automate this reduction as thoroughly as the factorization of entangled singularities. We 
will therefore discuss in detail the only case needed for Higgs hadroproduction, which is 
the l/si3 line singularity in the rapidity parameterization. We can avoid the I/S34 energy 
parameterization singularity, as we will show in the next section. We will discuss how to 
handle the topology I/S13/S24/S34, which is the simplest topology that can not be rotated 
away from a line singularity form. For notational simplicity, we will suppress any possible 
numerator and the measurement function Fj for this structure. 

We begin by combining the expressions for S13, S24, and S34 in Eq. ljU^ with the phase 
space in Eq. flU^ . We obtain 
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and Ai,A2 defined in Eq. flMj) . We note that the line singularity occurs when A4 = and 
Ai = A2. The first step is to separate these two requirements. To do so, we make use of 
the freedom to map singularities nonlinearly to the unit hypercube; so far, we have only 
used hnear transformations such as the one used for S13: S13 = — A4(si3 — Si^) — Si^. The 
nonlinear mapping is achieved by making the additional variable change 
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where we have relabeled A4 A4 for notational simplicity. The Ai = A2 line has been made 
manifest in the absolute value in the last line of this equation. As noted in Eq. ()(j5|l . we can 
describe this line by the parameter A2. If we split the A2 integration into the two regions 
[0, A2] and [A2, 1], we can force this singularity to always occur at a point on the boundary 
of phase space; this will reduce it to an entangled form, amenable to sector decomposition. 
We denote the integrals over these two regions by Iq and Ji, so that / = Jq + /i. In Jq we 
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perform the variable change A2 = A2/A2, while in Ii we use X'2 — (A2 — A2)/(l — A2). We 
obtain (after relabeling A2 ^ A2) 
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The terms with A^*^ and ^2*'' contain singularities entangled in Ai and A2 that can be 
extracted with sector decomposition. Although the above expressions are somewhat com- 
plicated, the simple, automated algorithm described above can separate and extract all 
singularities with minimal human intervention. We note that after splitting Iq and /i into 
the ranges [0,1/2] and [1/2,1] for Ai, A2, A3, and factorizing all singularities, we obtain 
twenty sectors for this topology, each with a different mapping between the Sab and the 
Aj. It is clearly advantageous to avoid line singularities by rotating them into entangled or 
factorized forms when possible. 



VIII. TOPOLOGIES FOR HIGGS HADRO-PRODUCTION 



We have discussed how to extract singularities from all terms that appear in the double 
real emission corrections to Higgs production at NNLO. We must now explain what singular 
structures appear in the matrix elements, what parameterizations we use to deal with them, 
and how the numerator structures arc handled. We consider the topologies that appear in 
9{Pi) + 9iP2) H{pfi) + giPs) + g{p4)', all other partonic channels contain only a subset of 
these topologies. 

There is a certain amount of ambiguity in how to map terms in the matrix elements to 
different topologies. For example, the matrix elements can contain two sets of terms with 
denominators 1/52/1/534 and 1/sih/ S23/S34. There arc two possible ways of deahng with these 
terms. We can treat them separately, or we can rotate pi ^ p2 in the first and combine 
them over a common denominator. Each approach offers potential advantages. Combining 
them allows us to extract singularities from only a single structure, 1/s 1/1/523/534, which 
reduces both the number of singular regions in Aj space, and the expression size. The first 
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term is singular as A3 — 1, while the second is singular as A3 ^ 0; the combined expression 
is singular only when A3 — >^ 0. However, combining these terms can complicate numerical 
cancellations between terms in the matrix elements. As noted, to combine these terms 
requires rotating pi ^ p2 in I/S2/1/S34, which changes the A3 — 1 limit to A3 0. Suppose 
there is a numerical cancellation between the I/S2/1/S34 term and another component that 
occurs locally in Aj-space when A3 ^ 1. After rotation, this cancellation will only occur 
globally after integrating over A3. 

The point of this example is to show the types of possibilities that exist when we map the 
matrix elements to topologies. It is not clear to us at the moment which way is "better": 
which leads to smaller expressions and to more numerically stable results, etc. In the dis- 
cussion that follows, we focus on the parameterization- independent features of the possible 
mappings. 

• We will note which topologies are necessarily of an entangled or a line-singular form. 

• We will indicate which topologies are "soft", i.e. singular in the z —>■ 1 limit. These 
only occur for the partonic channels with a gg initial state: gg ggH and gg —>■ qqH. 

• We will discuss the topologies with "quadratic" singularities; we will explain later in 
detail what this means. 

We will also try to explain what techniques we found useful in organizing our calculation, 
and will attempt to discuss the possible highlights and disadvantages of different methods. 

We first discuss how we deal with the numerator structures of each topology. Each term 
that appears in the matrix elements has the form 

Sah^cd • ■ • 

where the numerator M contains polynomials of the invariant masses in addition to mea- 
surement functions Fj with various arguments. The denominator characterizes the topology 
because it is insensitive to the exact particle content of the theory, and only depends on the 
global structure of the process under consideration. Identical topologies will appear for all 
2 — > 3 processes with one massive particle in the final state. On the contrary, numerators 
are process specific; for example, in a theory with only scalar particles, all the numerators 
are equal to one. It is therefore advantageous to treat the numerators and denominators of 
the matrix elements separately. 

One way to deal with the numerators in a numerical program is to explicitly substitute the 
expressions for the Sab in terms of the Aj for each topology, so that everything in our Fortran 
code is written in terms of A,. However, we found that it is better to define auxiliary 
functions ^^^(Aj) to keep N written in terms of the Sah in our numerical routines. We 
introduce the explicit formulae for the Sab only for the denominators, which is required in 
order to extract singularities. To show what appears in our numerical code, consider the 
topology A/'[sab(A2, . . ■)]/ 82^. This is singular as A2 — > 0, as is clear from Eq. (jUH|) . so we get 
terms such as N[sab{^2i ■ ■ •)]/[^2] + - What we write to our code is 

A/'k.(A2,...)]-ArM0,...)] ^ 
A2 

The terms Sab(A2, . . .) and ^^^(O, . . .) are obtained numerically by calls to a Fortran function. 
This drastically reduces expression sizes; for example, referring to Eq.(|U^. we see that if a 



25 



term has S13 in the numerator, it is much cleaner to keep it written as si3(Aj) rather than 
exphcitly introduce the lengthy functional form. Since all topologies depend upon a small 
set of Sab, optimization of the code is appreciable. As mentioned above, the denominator 
structures are universal; the same ones appear in Higgs hadroproduction, electroweak gauge 
boson production, bb ^ H production in supersymmetric theories, and a host of other 
phenomenologically interesting 2 — 1 processes. Structuring our numerical code by treating 
the denominators and numerators of the matrix elements differently allows us to study other 
processes simply. 

We are now ready to discuss the topologies that appear in Higgs hadroproduction. We 
first introduce some notation. We refer to a topology with four invariant masses as a 
"highest-level" topology, as this is the maximum number of scalar products that can occur 
in the denominator for the double real emission contributions. When we remove an invariant 
mass from a topology, we refer to this as a "sub-topology" of the original term. We remind 
the reader of the three classes of diagrams we introduced in the section on phase-space 
parameterizations. We will refer to them as Ci,C2,C3, and interferences between them as 
Ci X C2, etc. Only topologies independent under rotation will be discussed; for example, if 
we present I/S13/S24, we will not consider l/si4/s23- 

• The only highest-level topology that necessarily contains a line singularity is 
^/si3/s24,/s34^/sih. It occurs in C2 x C3. It is also a soft topology, and therefore 
begins at We found it best to map this term into the rapidity parameterization; 
with this choice, we obtain 20 sectors, as we noted in our discussion in the previous 
section on line singularities. If we had chosen instead the energy parameterization, 
more sectors would be required to factorize the singularities present in sih- This is the 
most difficult topology to evaluate numerically; the singular 2; — > 1 limit creates a large 
number of plus distributions at the finite level. It has a subtopology I/S13/S24/S34 that 
also has a line singularity; its remaining subtopologies do not. 

• There are several highest-level topologies that necessarily contain entangled singular- 
ities. The two most complicated of these are I/S13/S23/S1/1/S2/1 and l/si3/s24,/sih/s2h, 
which appear in Ci x C3. They are both soft, and begin at If we evaluate them 
using the energy parameterization, we can avoid line singularities associated with S13. 

• The other highest-level topology that is necessarily of entangled form is 
I/S13/S23/S24/S1/1 (another one that requires a special discussion will be considered 
later). It is not soft, and so only begins at It appears in Ci x C3. We can again 
avoid line singularities by using the energy parameterization. 

We must now discuss topologies with "quadratic" singularities. To explain this concept, 
we will discuss the simplest example where quadratic singularities occur. Consider the sub- 
topology l/s23/su/sih in the rapidity parameterization. When we combine the expressions 
in Eq.(jnSl) with the phase space in Eq. ljU^ . we find that this topology scales as Aj^"^"; 
the topology is quadratic in A3. Singularities this strong are unphysical. If we were to 
regulate them with a mass m rather than with e, we would find a rather than a 

In(m^) singularity, which can not occur. Furthermore, the plus distribution expansion in 
Eq. (|22j) is only valid for linear singularities. When we find quadratic singularities, we always 
have a term in the numerator of the topology that renders the behavior acceptable. In 
this example, we must find that the numerator scales as A3. However, finding this scaling 
requires introducing the explicit expressions for the invariant masses in terms of the A,. 
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We want to do as little of this as possible in order to keep term sizes small. Topologies 
with quadratic singularities in A3 are the simplest to handle. We map all such cases to the 
rapidity parameterization. Referring to Eq. ljUHj) . we observe that all the invariant masses 
have a simple, overall scaling in A3, together with a more complicated term. For example, 
Sih ~ A3, Sih ~ (1 — A3) Si3 ~ (1 — A3), etc. We find that introducing this overall factor is 
sufficient to cancel quadratic singularities; the terms in brackets that appear in these scalar 
products never enter the cancellation, and can be kept implicit. This allows us to keep the 
term sizes small. 

The remaining quadratic topologies, i.e., those with a quadratic singularity in a variable 
other then A3, require a slightly more involved procedure. There is typically a combination of 
invariant masses in the numerator for which we must introduce an explicit parameterization 
in order to regulate the singularity. 

• In the interference Ci x Ci, we find the topology l/sl^/s24^. This is factorizable in the 
energy parameterization; however, referring to Eq . (j45j) . we observe that it scales as 
(1 — A2)~^^^A3"^~^. The relevant numerator structure that regulates these quadratic 
singularities is 

(S14S23 - S34)' cx (1 - A2)A3(1 - z)\ (83) 
The (1 — factor that appears ensures that this topology is not soft. 

• Two other similar quadratic topologies that appear are l/s24/s^^ and 1/324/513/51^,. 
The first comes from C3 x C3. It is factorizable in the rapidity parameterization, where 
we find that it scales as (1 — A2)~^~^. The second topology comes from Ci x C3, and 
necessarily contains entangled singularities. We can avoid a line singularity in the 
energy parameterization, where we find that it also scales as (1 — A2)~^~^. They are 
both regulated by the same numerator structure as above, (S14S23 — 534)^. Neither 
topology is soft. 

• The final two quadratic topologies appear in C2 x C2: ^/sljslf^ and 1/ sl^^/ sih/ S2h- 
They are factorizable in the rapidity parameterization, where we find that they scale 
as Af^"*^. The required numerator structure is 

{ZS23 - SisSih - S23Sih - S23Y oc AiA3(l - A3)^(l - z)^, (84) 

which also ensures that the topology is not soft and that quadratic singularities in Ai 
cancel. 

All of the remaining topologies that are needed for Higgs hadropro duct ion can be mapped 
directly to a factorized form in either the energy or rapidity parameterization. We find the 
following highest-level topologies: 

• ^/ S23 /s24,/slf^ in C3 X C3, which we map to the rapidity parameterization; 

• l/'^23/s34/si/j in C2 X C3, which we map to the rapidity parameterization; 

• ^/ S23/ S3i/ sih/ S2h in C2 X C3, which we map to the rapidity parameterization; 

• I/S13/S14/S23/S24 in Ci X Ci, which we map to the energy parameterization. 

Only the final two topologies in this list are soft. We can directly apply the expansion in 
Eq.(I221) to these terms to extract the phase-space singularities. 



27 



IX. RESULTS 



In this Section we describe phenomenological results obtained using our calculation. Our 
primary goal is to illustrate the range of observables that can be studied using the approach 
discussed in this paper, not to perform a comprehensive study of the Higgs boson signal at 
the LHC A state-of-the-art phenomenological analysis of the Higgs boson signal is required 
for two reasons. First, it can be used to optimize cuts, which is important for enhancing the 
signal-to-background ratio. Second, it is useful for exposing uncertainties in the theoretical 
prediction for the Higgs boson signal, including truncation of the perturbative expansion at 
NNLO, choices of the factorization and the renormalization scales, and imprecise knowledge 
of parton distributions and as{Mz). Although such a study is beyond the scope of this 
paper, we hope that the examples presented below are sufficiently convincing to demonstrate 
that our approach permits a computation of arbitrary observables related to the reaction 
pp H + X Higgs decay products In particular, we present below the fully realistic 
NNLO cross sections for the Higgs boson signal in the di-photon decay channel, where the two 
photons in the final state satisfy all the selection criteria (cuts on photon pseudorapidities, 
transverse momenta, and geometric isolation from significant hadronic activity) used by the 
ATLAS and CMS collaborations. 

We remind the reader that we work in the limit of an infinitely heavy top quark, in which 
the Higgs boson coupling to gluons is point-like. However, all the results presented in this 
Section are rescaled by a factor 

^{mt) = J r- (85) 

(^Loim = oo) 

This rescaling accounts for the effects of the finite top mass exactly for LO cross sections, 
and provides an approximate description of the top mass dependence at higher orders. For 
NLO calculations, this is known to be an excellent approximation, and we expect it to work 
well at NNLO also. 

All the results that we present in this Section can be divided into two categories, de- 
pending on whether or not the decay of the Higgs boson into two photons is considered. 
We first study the simpler case, in which the Higgs boson is treated as a stable particle, 
and investigate the various kinematic distributions of the Higgs. We then turn to the anal- 
ysis of the Higgs signal in the di-photon channel, and present the NNLO cross sections and 
the distributions in the pseudorapidity difference \ri'^'^ — rp''^\/2 and the average transverse 
momentum (p]^'^ -|-p]^'^)/2 of the two photons. 

The rapidity distribution of the Higgs boson with mass = 120 GeV is shown in Fig. El 
at LO, NLO, and NNLO in QCD perturbation theory. The bands are obtained by equating 
the factorization and the renormalization scales and then varying /x in the range 

mh/ 2 < yu < 2mh. We use the parton distribution functions (pdfs) from the MRST2001 
set j7C| . and employ mode 1 as the default mode for pdf evolution. Higher order QCD 
corrections to the Higgs rapidity distribution exhibit behavior similar to the corrections to 
the inclusive cross section; the large increase in dcr/dY from LO to NLO, is not followed by 
a similar large increase from NLO to NNLO. The NNLO corrections are uniform over the 
central rapidity interval |y| < 2, and can therefore be obtained to good approximation by 
re-scaling the NLO rapidity distribution by a universal, rapidity-independent factor. 

For a Higgs boson heavier than ~ 140 GeV, the decay H — > W^W~ becomes the dom- 
inant decay mechanism, making the Higgs branching ratio into two photons quite small; 
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FIG. 3: Bin- integrated Higgs boson rapidity distribution at the LHC. The bands indicate the scale 
choice m/i/2 < /i < 2mh- 

consequently, the two photon signal can not be used as the primary trigger for the Higgs. 
Searching for the Higgs boson in the W~^W~ decay mode requires the introduction of addi- 
tional cuts to suppress the background due to the production and subsequent decay of a pair 
of top quarks. Since the hadronic jets in top pair production have, on average, larger trans- 
verse momenta than hadronic jets in Higgs hadroproduction, the significance of the Higgs 
signal can be enhanced by imposing a jet veto on the recoiling hadronic system |7ll. It^. In 
Fig. m we present the rapidity distribution of the Higgs boson with mass rrih = 150 GeV, 
when all jets in the final state of the reaction pp —>■ H + X are required to have transverse 
momenta smaller than Px^veto — 40 GeV. The jets are identified with the cone algorithm, 
using a cone size R = OA. 
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FIG. 4: Bin-integrated Higgs boson rapidity distributions at the LHC with a jet veto applied. 
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It is instructive to compare both the magnitude of the perturbative corrections and the 
residual scale dependence of the Higgs rapidity distributions with and without the jet veto. 
Although Figs. I3I4I refer to different masses of the Higgs boson, a qualitative comparison 
is still possible; as is known from the NNLO computation of the inclusive cross section, 
the relative magnitude of the NNLO corrections does not depend strongly on the Higgs 
boson mass. From Figs. I3I4I we observe that the application of the jet veto leads to a 
better perturbative stability of the rapidity distribution; the perturbative corrections are 
smaller when the jet veto is applied, and there is a complete overlap of the NNLO scale 
dependence band with the NLO band, for most rapidities. This comparison shows that large 
and perturbatively unstable contributions to the Higgs boson hadroproduction cross section 
are related to kinematic configurations where the Higgs is produced with large transverse 
momentum. This feature does not imply a breakdown of the perturbative expansion; it 
merely reflects the fact that for Higgs hadroproduction, the leading order partonic process 
is gg —>■ H, so that at LO the Higgs boson is produced with zero transverse momentum. 
Therefore, if Higgs boson production with p^^^^ 7^ is considered, our NNLO computation 
includes just two terms in the perturbative expansion in the strong coupling constant, and 
is therefore a NLO computation. 

The excellent convergence observed for the jet- vetoed rapidity distribution with Px veto = 
40 GeV is partly accidental. This can be seen from Fig. where we show the dependence 
of the Higgs boson production cross section on the value of a veto on the transverse energy 
of the Higgs, PT,veto- While the LO cross section obviously does not depend on PT,veto, both 
the NLO and the NNLO cross sections exhibit a significant dependence on this parameter. 
It is interesting to observe that the NLO cross section reaches its asymptotic value much 
faster than the NNLO one; this is related to the fact that the average of the Higgs boson 
increases from NLO to NNLO. For example, for rrih = 150 GeV, the average transverse 
momenta of the Higgs boson at NLO and NNLO are = 37.5 GeV and (p^^'"^) = 

44.6 GeV. It is clear from Fig. El that the choice PT,veto = 40 GeV minimizes the NNLO 
QCD corrections for nih = 150 GeV. However, other choices of pT,veto also lead to only 
small differences between the NLO and the NNLO cross sections, indicating an improved 
convergence of the perturbative expansion for all veto choices. 
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FIG. 5: The Higgs boson production cross section at the LHC as a function of px.veto- 
Another interesting feature of the NNLO result is the dependence of the jet- vetoed cross 
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section on the cone size R. At NLO, every jet is associated with a massless gluon or quark, 
so there is no cone-size dependence of the cross section. At NNLO, when two partons are 
close in rapidity and azimuthal angle, they are combined into a single jet. Thus, a cone- 
dependence of the observable cross section appears. As can be seen from Fig. IHl the cross 
section decreases when the cone size is increased over a large range of i?; however, at large 
values of R the cross section starts to increase again. This is a consequence of the fact 
that, originally, combining two energetic partons into a single jet increases its transverse 
momentum. However, after the cone size becomes so large that two partons which are back- 
to-back in the transverse plane are combined to form a single jet, the momentum of such 
a jet becomes smaller than the momenta of the individual partons. This effect drives an 
increase in the cross section for R > 2.5. 
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FIG. 6: The Higgs boson production cross section at the LHC with a jet-veto of pT.veto = 40 GeV, 
as a function of the cone size R. 



Since the rapidity of the Higgs boson is not probed, the results of Fig. El ca n also be 
obtained fn\ by subtracting the NLO contributions for H + 1 iet production jSOl |6^ from 
the NNLO total cross section. We have compared our results with those in Refs. [S^i, ^flj , 
and have found excellent agreement. 

We now discuss the Higgs boson signal in pp ^ H + X — 77 -|- X. An interesting 
observable is the total di-photon production cross section, subject to a number of cuts on the 
kinematic variables of the two photons. These cuts, designed to enhance the signal over the 
prompt photon background, include restrictions on the photon transverse momenta and the 
rapidities, as well as isolation cuts. In particular, the photons are required to have transverse 
momenta p^^^ > 40 GeV and p^f^ > 25 GeV; they must also be produced in the central 
rapidity region |?7| < 2.5 i^l- The isolation cuts are designed to reduce contamination of the 



signal by the poorly known fragmentation contributions. We require that a photon candidate 
does not have an additional transverse energy which is greater than -Ex.min = 15 GeV 
deposited within a cone around it of radius R^g = \J {v ~ ViY + (0 ~ 'P-yY = We will 
refer to this combination of cuts as the "standard photon cuts" in what follows. 

In Table H] we present the cross section for pp —>■ H + X — > 77-l-X, divided by the 
branching fraction of H —y 77, with the standard cuts imposed on the photons. We give 
results through LO, NLO and NNLO in perturbation theory, for different Higgs boson masses 
and for two choices of the renormalization and factorization scales. The NNLO results are 
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TABLE I: The cross section for pp-^H + X—^jj + X, divided by the branching ratio Br{H — > 
77), at ^/s = 14 TeV with the standard cuts apphed to the photons, for different values of the 
Higgs mass. 
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TABLE II: Comparisons between the cut and inclusive cross sections for different Higgs masses. 
The second column contains the ratio of the NNLO cross section with the standard cuts over the 
inclusive cross section, while the third column contains the ratio of cut and inclusive results for 
the ivT-factor K^"^^ = (Tnnlo/<7nlo • We have set n = mh/2. 

accurate to 1%, while the LO and NLO numbers are accurate to 0.1%. The perturbative 
corrections to the di-photon signal follow the pattern of the corrections to the inclusive 
Higgs production cross section. We note the apparent convergence of the di-photon signal, 
with the NNLO corrections increasing the NLO result by up to ~ 20%, depending on the 
value of n and the Higgs boson mass. The NNLO and NLO scale dependence bands overlap 
significantly, and the NNLO scale dependence is reduced by a factor of two compared to the 
NLO dependence. 

It is interesting that the NNLO corrections are much smaller when the scale pL — mh/2 is 
chosen. This is in accord with the observation that the typical transverse momentum of the 

Higgs boson is parametrically smaller than the Higgs mass m^; because of the rapid increase 
of the gluon pdf at low x, the Higgs bosons are produced close to threshold. Since the 
factorization scale should be chosen close to a typical transverse momentum of the Higgs, 
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selecting /i = mh/2 expedites the convergence of the perturbative expansion. Fig. [7| shows 
the dependence of the Higgs signal on the choice of the scale for nih = 120 GeV. As expected, 
we find that the NNLO perturbative corrections are small for /i ~ 40 — 50 GeV. We note 
that the threshold resummed results for the Higgs hadropro duct ion cross section js^ agree 
very well with the fixed order results for smaller scale choices such as /i ~ Tnh/2, while 
they differ from the fixed order results by up to several percent for larger scale choices. It 
appears, from the stability of the perturbative series and the agreement with the resummed 
result, that /x ~ mh/2 is a better scale choice for Higgs hadroproduction. 

pp^77+X 
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FIG. 7: The Higgs channel pp — > H-\-X 77 + X with the standard cuts imposed on the photons, 
as a function of scale choice. 

In Table im we present several comparisons between the inclusive NNLO cross section, and 
the NNLO cross section computed with the standard cuts. Our goal in doing this is two-fold: 
to understand how the cross section is affected by the standard cuts, and to see how well the 
realistic cut cross section can be approximated if only the inclusive NNLO result is known. 
We adopt = mh/2 for these comparisons. These results are valid to approximately 2%. In 
the second column of Table ITTl we have presented the ratio of the NNLO cross section with 
the standard cuts over the inclusive NNLO result. The reduction of the cross section caused 
by the cuts ranges from 40% for mh < 120 GeV to 30% for mh > 150 GeV. We note that 
most of this reduction comes from the p± and r] cuts; the isolation cuts decrease the ratio 
by less than 3%. This is expected; there is no cross section enhancement when a parton 
is emitted along the photon direction, so this phase-space region contributes minimally to 
the total result. We note that the cuts become less effective at larger mh, i.e., the ratio 
increases. For larger Higgs masses, the average photon p± increases, and therefore more 
events pass the cuts. 

In the third column of Table ITTl we present the ratio of the K-factor K^'^^ = ctnnlo/'^nlo- 
This is interesting for the following reason. Suppose only the differential NLO cross sec- 
tion and the inclusive NNLO result are known. The best approximation for the exact 
NNLO differential result would then be da^^^Q = da^i^o x K^^l, where K^^l is defined with 
the inclusive cross sections. Calculating the cut cross section with this distribution gives 
'^NNLo''^'^* — ^NLO ^ ^inc- The ratio of this result with the exact NNLO cross section with 
the standard cuts imposed is c''nnlo /'^nnlc) '^"^ — ^aL/^inc'i the deviation of this ratio from 
unity measures the error made by using da^^]/^Q to approximate the actual differential cross 
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section at NNLO. We see that this deviation is less than about 5%. 

In order to optimize the experimental cuts, it is desirable to have a good understanding 
of the kinematic distributions of the photons, since they can provide good discriminators 
between the signal and the background. While we do not discuss cut optimization in this 
paper, we present two differential distributions that illustrate the range of observables that 
can be studied using our calculation. In Fig. |H1 the p± = {p]^^ +p']^^)/2 distribution is 
shown for rrih = 120 GeV. We observe large perturbative corrections close to the kinematic 
boundary at leading order, p± < mh/2 = 60 GeV, where resummation of large logarithms 
is required. However, the presence of a large peak near the LO kinematic boundary ap- 
pears to be a reliable result, as it appears without drastic modification at both NLO and 
NNLO. Since the background should not contain any such feature, this is potentially a useful 
discriminating variable. 
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FIG. 8: The p± = {p'Y^ +pj^'^)/2 distribution for the di-photon Higgs signal at the LHC. 



In Fig. ini we present the distribution of the pseudorapidity difference Kj = jr^^'^ —ri'^''^\/2 
between the two photons. This distribution is interesting since a similar distribution from the 
prompt photon production background is flatter; this information can be used to enhance 
the statistical significance of the Higgs signal |4Q]. From Fig. IHl we see that the peak at 
|j^7.i _ j^7.2| = is also present when the NNLO effects are included. 

The results presented in this Section are for the di-photon Higgs signal. There are other 
Higgs decay modes that are of significant interest. In particular, for moderately heavy Higgs 
bosons, decays into ZZ —>■ 41 and W~^W~ — >■ might provide suitable channels for 

discovery. Since our calculation retains all the information about the Higgs boson kinematics, 
it is in principle straightforward to include Higgs decays into arbitrary final states. However, 
in reality, some care must be exercised to generate the final-state decay efficiently, especially 
for decays with high multiplicities and sophisticated cuts. We plan on adding additional 
Higgs decay channels to our code in the future. 



X. DESCRIPTION OF THE FORTRAN CODE 

In this Section we describe a FORTRAN program, FEHiP, which we have written to obtain 
the results described in the previous Section. As can be seen from the examples presented 



34 




0.00 0,25 0.50 0.75 1.00 1.25 1.50 



Ys 

FIG. 9: The distribution in the pesudorapidity difference of the two photons, Yg = \r]J — r/2 |/2. 

there, FEHiP computes the cross section for Higgs boson production in hadronic coUisions 
through NNLO in perturbative QCD in the presence of an arbitrary measurement function. 
The decays of the Higgs are treated in the narrow width approximation. At present, only the 
Higgs decay into two photons is included; other decays can be added. This Section provides 
instructions for using the code. 

A. Download and Compile 

The code can be downloaded from j?^. To compile the code, first uncompress the file 
fehip.tar.gz and run the "make" script: 

\home\user> tar -zxvf fehip.tar.gz 
\home\user> cd FEHIP 
\home\user\FEHIP> make 

We have sussesfully compiled the code on Linux systems using the GNU g77 compiler. To 
compile the code on other platforms the user should modify the file makefile. To run the 
code, execute the program f ehip: 

\home\user\FEHIP> fehip 

The program performs multidimensional numerical integrations using the adaptive Monte 
Carlo algorithm Vegas j?^; we use its recent implementation in the CUBA library [t^. The 
current version of this library, Cuba-1.0, is included in the distribution of FEHiP, and 
compiles automatically with the above make script. For future updates of the Cuba library, 
we refer the user to [75|; the directory FEHiP\Cuba-l . 0\ and the file makefile should be 
updated appropriately. 

B. Basic parameters and Usage 

The basic parameters that are used by FEHiP include the mass of the Higgs boson, the 
type and energy of the hadron collider, the factorization and renormalization scales, and 
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the order in perturbation theory through which the result is to be computed. The code 
also requires the mass of the top quark and the value of the Fermi constant to compute 
the normalization of the Higgs boson production cross section (cf. Eq. ()85|) ). This input is 
provided by the user in the file input.txt. 

A prototype input file has the following format: 

'Mass of the Higgs boson (GeV) = ' 120d0 



'Collider (pp=0, ppbar=l) = ' 

'CMS collision energy (GeV) = ' 14000dO 



'Factorization scale (GeV) = ' 120d0 
'Renormalization scale (GeV) = ' 120d0 



'Top-quark mass (GeV) = ' 175d0 

'Fermi constant = = ' . 0000116639d0 



'Final-state (O=no-decay, l=diphoton) = ' 1 
'Branching ratio = ' IdO 



Perturb. Order (0=L0, 1=NL0, 2=NNL0) = ' 1 



'Output File =' 'result.dat' 

The user should provide values for the following variables: 

• The mass of the Higgs boson. This is a double precision variable that sets the 
value of the Higgs boson mass nih, in GeV. 

• Collider. An integer variable that defines the type of collider. For proton-proton 
collisions, the value must be set to 0; for proton-antiproton collisions, the value must 
be set to 1. 

• Energy. A double precision variable for the energy of the collider, in GeV. 

• Factorization scale. A double precision variable for the value of the factorization 
scale fip, in GeV. 

• Renormalization scale. A double precision variable for the value of the renormal- 
ization scale ^R, in GeV. 

• Top-quark mass. A double precision variable for the value of the top-quark mass 
mt, in GeV. 

• Fermi-constant. A double precision variable for the value of the Fermi constant Gp, 
in GeV^l 

• Final-state. An integer variable that determines the decay mode of the Higgs boson. 
The value corresponds to the "no-decay" option. In this case it is assumed that 
the Higgs boson four momentum is fully reconstructed and is therefore an observable. 
With the "no-decay" option, it is possible to impose cuts on the transverse momentum 
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and the rapidity of the Higgs boson, and to study jet-clustering and jet-veto cuts for 
the additional radiation (of up to 2 partons) in the final-state. However, no information 
about the kinematics of the Higgs decay product is provided. 

The value 1 corresponds to the "diphoton decay mode" . In this case, the Higgs decays 
into two photons in the narrow width approximation, and the kinematics of the two 
photons in the final state is fully reconstructed. Hence, arbitrary cuts on the photon 
momenta can be imposed. 

We note that the results for the "no-decay" mode can always be obtained from the 
"decay" mode. Indeed, if no restrictions on the momenta of the individual photons 
are imposed, the two options should provide identical results. However, if one is not 
interested in the kinematics of the Higgs decay product, it is beneficial to declare 
the "no-decay" option explicitly, since the dimensionality of the numerical integration 
performed is then lowered. 

Finally, we note that in the future, additional decays of the Higgs boson will be added 
to the program. 

• Branching ratio. A double precision variable for the value of the branching fraction 
of the Higgs boson decay into the selected final state. The program computes the 
Higgs production cross section in the narrow-width approximation and multiplies the 
output by the 'Branching ratio'. The value of the branching fraction for a particular 
decay of the Higgs with a certain mass should be supplied by the user as an input. 
Theoretical predictions for the branching fractions of the Higgs can be obtained with 
the program H DECAY [t^. 

• Perturbative Order. An integer variable which sets the order through which the 
perturbative expansion of the Higgs boson production cross section is computed. The 
values 0, 1, and 2 correspond to the LO, NLO, and NNLO cross sections, respectively. 

• Output File. A character string variable for the file name where the results of the 
calculation are written. The output file has the following format: 



proton-proton collider 

CMS collision Energy = 14000. 



Higgs boson mass = 


120. 


Factorization scale 


120. 


Renormalization scale = 120. 


Top-quark mass = 


175. 


Fermi constant = 


1.16639E-05 


H — > gamma gamma 




Branching Ratio = 


1. 


NLO Cross Section 




Strong coupling = 


0.114230586 
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=========== RESULT 

Sigma = 38.4079658 +/- 
chi square probability = 



0.0362102225 
2.46758092E-05 



It contains a description of the input parameters and the result (sigma) for the cross 
section in picobarns with the corresponding statistical error from the Vegas integration. 
It also provides the x^— probability (prob variable from the Vegas routine in CUBA), 
which is an indicator of the reliability of the Monte-Carlo integration. 

A number of other necessary program settings and options are not included in the user 
controlled input.txt file. Some important settings can be found in the main program file 
higgs .F; we have set them to values that are reasonable for practical purposes, according to 
our own experience with the code. In this file, we define the parameters for Vegas, compute 
the normalization of the cross section, and invoke the integration routine. For example, we 
require a precision (the epsrel variable) of 1% at NNLO and of 0.1% at lower orders with 
an absolute error above [10~^ x Br] ph. We note that FEHiP also writes the iteration- by- 
iteration results produced by Vegas to screen, so the user can track the output; this can 
easily be redirected to a separate file. This can be suppressed by modifying the appropriate 
fiag in higgs. F; we refer the reader to the CUBA documentation for a discussion of the 
various output options. 

FEHiP uses the MRST 2001 LO, NLO, and NNLO parton distribution functions [tJ. 
These are provided with the distribution of FEHiP. We note that it is extremely inefficient 
to directly call the pdf routines. Since the factorization scale and the range of the Bjorken 
variable x are fixed during each run of the program, we only need a limited amount of 
information about the pdfs. To optimize the code, we adopt the following strategy. When 
the code is initialized, FEHiP first calls the pdf routines for a given value of the factorization 
scale and computes the pdfs for 1000 different choices of x. The distribution of these points 
is not uniform; we increase their density for x values where the pdfs change rapidly. We then 
use interpolation with quadratic polynomials to connect the adjacent points. Consequently, 
each pdf is described by an array of coefficients of these quadratic polynomials which is kept 
in memory; this speeds up the execution of FEHiP enormously. We have tested that the 
approximate pdf values are accurate to better than 0.05% for all x and hf values relevant 
for Higgs production. 

To use a different set of pdfs, the user should replace the mrstpdfs.F file with the 
appropriate set, and modify the calls to the mrstlo, mrst2001, and mrstnnlo subroutines 
in the file f itpdf .F. In addition, the user should change the value of as (Mz), which is used 
as an initial condition for the LO ,NLO, and NNLO evolution of the strong coupling, to the 
value which is consistent with the fitting of the new pdf set. The variables asZlo, asZnlo, 
and asZnnlo in the file higgs. F correspond to the fitted values for the MRST 2001 LO, 
NLO, and NNLO pdfs. 

C. Experimental cuts 

The most complicated user input is the definition of the experimental observable; this 
requires imposing cuts on the phase-space of the final state. We describe here a double 
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precision function, named constraint, in the constraint .F file, which can be used for this 
purpose. 

The constraint function corresponds to the observable function Fj for the final state. The 
NNLO parton multiplicity is assumed: 

partoni(pi) + parton2(p2) H{ph) + parton3(p3) + parton4(p4). (86) 

All final-state configurations with lower multiplicities have been mapped to this "maximal" 
final state by introducing additional partons with zero momentum when needed; this allows 
a uniform introduction of cuts for both real and virtual corrections. 
The routine uses the following (local and global) arguments : 

• the Vegas-generated variables: var; 

• the dimensionality ndim of the Vegas integration; 

• the squared mass of the Higgs boson, m2 = m^, and the CMS square energy of the 
colliding hadrons scm = s; 

• the Biorken variables xl, x2 and the ratio z = ; 

•> ' SX1X2 ' 

• the variables sl3,s23,sl4,s24 with Sif — — i = 1,2, / = 3,4, and the variable 

• the variables slv,s2v with Siy = ^ = 1,2; 

• the variables s3v,s4v with s„/ = ^^7^7^, / = 3, 4. 

From these variables we can reconstruct the components of the momenta of all particles in 
the final state, and impose cuts which define the observable. We can also generate (using 
the var array) the required additional phase-space variables for the subsequent decay of the 
Higgs boson. 

After the appropriate variable declarations at the beginning of the routine, we proceed 
by generating phase-space variables for the two photons coming from the decay 

H{ph)^-fa{qa)+%{qb), (87) 

using a convenient parameterization in terms of invariant masses formed from the photon 
momenta Qa^b and the initial state partonic momenta pi^2- The routine, for each Vegas event, 
computes the following quantities. 

• The transverse energy of: (1) the Higgs boson, pT; (2) parton 3, ptS; (3) parton 4, 
pt4; (4) 7„, ptga; (5) 7^, ptgb. 

• The energy of: (1) the Higgs boson. En; (2) parton 3, EnS; (3) parton 4, En4; (4) 7^, 
Ena; (5) 7;,, Enb. 

• The momentum along the beam axis of: (1) the Higgs boson, pZ; (2) parton 3, 
pz3; (3) parton 4, pz4; (4) 7„, pza; (5) 7b, pzb. 
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The above variables can be defined for all events, including the ones where unresolved 
(soft/collinear) partons and photons are generated in the final state. However, the pseu- 
dorapidity and relative angles of massless particles can be sensibly defined only if they are 
resolved. To discriminate, we introduce a cutoff for the transverse energy of photons and par- 
tons in the final state, ptbuf = O.OlGeV. We then use two flags, topf lag and photonf lag, 
to characterize the event according to how many partons and photons have a transverse 
energy above the cutoff. For every event we define: 

• topf lag = OdO: two unresolved partons in the final state, pt3,pt4 < ptbuf. 

• topf lag = IdO: one resolved and one unresolved parton in the final state, pt3 < 
ptbuf, pt4 > ptbuf or pt4 < ptbuf, pt3 > ptbuf. For such an event we compute 
the additional variables: 

— the transverse energy of the resolved parton: ptr. 

— the rapidity of the resolved parton: etar. 

• topf lag = 2d0: two resolved partons pt3,pt4 > ptbuf, and the Higgs boson with 
non-zero transverse energy, pT > 0. For this event we also compute: 

— the pseudo-rapidity of parton 3, etaS, and parton 4, eta4. 

— the azimuthal angle between parton 3 and parton 4: phi34. 

• topf lag = 3d0: corresponds to the (rare) configuration with two resolved partons, 
pt3,pt4 > ptbuf, and the Higgs boson with zero transverse energy, pT = 0. For this 
event we also compute: 

— the pseudo-rapidity of parton 3, eta3, and parton 4, eta4. 
We also define: 

• photonf lag = OdO: at least one unresolved photon in the final state, ptga < ptbuf 
or ptgb < ptbuf. 

• photonf lag = IdO: corresponds to two resolved photons ptga, ptgb > ptbuf, and 
the Higgs boson with non-zero transverse energy pT > 0. For this event we also 
compute: 

— the pseudo-rapidity of ■ja, etaa, and jb, etab. 

• photonf lag = 2d0: two resolved photons pt3,pt4 > ptbuf, and the Higgs boson 
with zero transverse energy, pT = 0. For this event we also compute: 

— the pseudo-rapidity of 7a, etaa, and 7b, etab. 

In order to study the isolation of photons, the routine computes the relative azimuthal angles 
of resolved photons with resolved partons. This can be done for the following combinations 
of flag values. 

• photonf lag = Id0,2d0 and topf lag = 2d0: The program computes the azimuthal 
angles for the pairs 
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— 7o, partong: angle phi3a. 

— 7b, partong: angle phiSb. 

— 7(j, parton4: angle phi4a. 

— 7b, parton4: angle phi4b. 

• photonf lag = IdO and topf lag = IdO: The program computes the azimuthal angles 
for the pairs 

— 7a, resolved parton: angle phira. 

— 7b, resolved parton: angle phirb. 

The above variables completely describe the phase-space for Higgs production and decay 
into photons through NNLO, and can be used to implement standard experimental cuts. 
We now describe the structure of constraint . F, and indicate how the reader modify the 
file to include other cuts. 

The section of constraint.? marked as STEP 1 performs the generation of the photon 
phase-space. STEP 2 computes the variables we have just described. Wc have placed 
the cuts required for the numerical results of this paper in the section labeled STEP 3: 
IMPLEMENTATION OF CUTS. These are: 

• the Higgs boson rapidity cut; 

• the Higgs boson px cut; 

• the jet- veto; 

• cuts on the pseudorapiditiy and pt of each photon; 

• an isolation veto on each photon; 

• a cut on the average pt of the photons; 

• a cut on the pseudorapidity difference of the photons. 

An event is accepted if the constraint function returns the value IdO, and rejected when 
it returns OdO. Before any cut is applied, the value of the function is set to IdO. Cuts are 
applied successively; each cut could potentially reject the event by setting the function value 
to OdO. To control which cuts should be apphed to the event, we use two flags: active, 
inactive. The user can choose combinations of the above cuts by modifying the appropriate 
flags, as described in constraint .F. 

The cuts that are programmed in STEP 3 of the routine serve as guiding examples for 
the user. When programming a new cut, some care is required for constraints on variables 
which are only defined when a final state particle is resolved; the user should provide that the 
constraint function returns an appropriate value for the events where the probed variable 
is not defined. In the section labeled STEP 4 of the constraint. F file, we have placed a 
"general" cut to guide the reader through this process. By modifying this section, it is 
possible to impose any desired constraint on the final state. 
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D. Other information 

To summarize, the user input for FEHiP is a set of parameters in the input.txt file, 
and the constraint function in the constraint . F file. To run the program, the user should 
compile using the make script and execute fehip. The make script should be executed every 
time that the constraint . F file is modified. 

We now describe some standard running times for the observables considered in this 
paper. We measure run times using the number of Vegas evaluations required to reach a 
given precision; on a 3 GHz PC, 10^ evaluations takes approximately 1 hour. The required 
number of evaluations for observables in the "no-decay" mode is typically fairly small. To 
reach 1% precision on the fully inclusive NNLO cross section at the LHC, about 1.5 — 2 x 10^ 
evaluations are needed. To reach 1% precision on a jet-vetoed cross section, about 8 x 10^ 
evaluations are needed. These numbers increase dramatically when the "decay" mode is 
used, and the standard cuts on the photons are imposed. For 1% precision on the pp —>■ 
if + X — i> 77 + X cross section with px, r], and isolation cuts imposed, 7 x 10^ evaluations 
are required. Computing bins for photon distributions to 2% or better precision usually 
takes over 10^ evaluations. 

Clearly, some observables require long running times. If the program is interrupted, the 
user can restart it by executing fehip. The CUBA library saves in a file the data of every 
Vegas iteration; this is used to restart the program from the last completed iteration. To 
force the program to restart from the first Vegas iteration, the user should type make clean 
and then execute fehip. 



XI. CONCLUSIONS 

We have discussed the calculation of the Higgs boson production cross section in hadronic 
collisions through NNLO in perturbative QCD. The kinematics of both the Higgs and the 
QCD radiation is kept exactly. This allows us to consider any decay of the Higgs, and impose 
arbitrary cuts on the final state. In this paper, we have focused on Higgs decays into two 
photons. We stress that this calculation provides the first example of a fully realistic NNLO 
calculation of the Higgs signal in the di-photon channel, where the two photons in the final 
state satisfy all the selection criteria (cuts on photon pseudorapidities, transverse momenta, 
and geometric isolation from significant hadronic activity) used by the ATLAS and CMS 
collaborations. 

To perform this calculation, we used the method of handling double real radiation de- 
scribed by us in js^, 3^|. We have further developed the approach in this paper, and have 



described its application to Higgs hadroproduction in detail. 

We have given a description of the FORTRAN code FEHiP which we used to obtain all 
the numerical results presented in this paper. The code can be obtained from [t^]. We 
hope that FEHiP will be used when state-of-the-art knowledge of the Higgs signal at the 
LHC is required. In particular, since arbitrary cuts on the two photons and accompaning 
hadronic radiation are allowed, the code can be used to optimize cuts to enhance the signal- 
to-background ratio. 

It was possible to obtain the results reported in this paper by further developing the 
approach suggested in [s^, 3^. Since this computation is far more complex than those 



considered in [3j, |36[, we believe that this is an important milestone for our approach. 
However, further development of the method is desirable. We give below a discussion from 
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this perspective. We indicate areas where more work is required, and describe possible 
sohitions to the problems of the method. We begin by listing attractive features of this 
approach. 

• Given a convenient parameterization of the phase-space for the real emission contri- 
bution to a given process, the singularities are extracted in an automated fashion; the 
human intervention required to achieve the final result is minimal. No classification 
of the various singular regions is required. Similarly, the cancellation of singularities 
is performed completely numerically, without the need for any analytic integrations. 
In principle, this method therefore provides an algorithm for the extraction and can- 
cellation of singularities through the N°LO order in perturbation theory. 

• Not only "real" singularities, but also some integrable ones are written in a factorized 
form. Consider the integral 



Although the singularity at the point a: ~ y ~ in /i is integrable, we can apply the 
algorithm described in Section VII to it anyway. Ii then becomes 



It is clear that the integrals I^' are now written in forms very convenient for numerical 
evaluation. In addition to extracting real singularities, our technqiue also smooths 
integrable singularities, thereby improving the numerical behavior of the integrand. 

• The complete kinematic information of the process is preserved. In principle, this 
enables us to use this result to construct a partonic level Monte-Carlo event generator. 

In spite of these attractive features, there are also some problems with this approach. 

• The most important drawback of the method is that it produces large expressions, 
which require long run-times for numerical evaluation and lead to difficulties with 
optimizing the code. Unfortunately, this is a natural feature of the approach; as 
we pointed out earlier, each time an entangled or line singularity is extracted, the 
expression size increases. However, it is useful to think of ways to ameliorate this 
behavior. It is instructive to recall how a somewhat similar problem was solved for 
NLO calculations. It was found that squaring matrix elements for unpolarized initial 
and final states is not a very practical option, since huge expressions are produced. 
Very compact expressions for scattering amplitudes can be obtained by working in the 
helicity basis. A similar approach can be attempted here. However, a given helicity 
amplitude contains terms that belong to different topologies, in the language of this 
paper. Since the separation into topologies is a crucial element of our approach, it is 
not clear how to to use helicity amplitudes efficiently. 

• By default, this approach might give more information than desired. All possible 
contributions to the cross section are automatically obtained. For example, in the 
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case of Higgs hadropro duct ion, we calculate H + 2 jets at LO, if + 1 jet at NLO and 
H + jets at NNLO. Imagine now that we want just the cross section of H + 1 jet 
at NLO. Obviously, the number of entangled and line singularities we must deal with 
should be smaller. However, a blind application of the algorithm in Section VII will 
automatically handle all singularities, including those needed for H + jets at NNLO. 
It should be possible to implement additional criteria in the algorithm to prevent 
this from happening; for example, a check on the Higgs pt could be included, or a 
counting of the number of lets could be imposed. In the calculation of e~^e~ — > 2 
jets at NNLO described in S^, 3^, the e+e^ 3 jets cross section at NLO can be 
obtained efficiently by including a call to the jet algorithm in the routine that handles 
singularities. However, a better understanding of these restrictions is needed before 
the approach can be used efficiently for processes with more complicated final states. 

Making a partonic level event generator out of our result is tedious, although possible. 
Because each sector corresponds to a different mapping of the invariant masses into 
A-variables, each sector must be generated as a separate channel in a multi-channel 
Monte-Carlo. Our result for the Higgs hadroproduction cross section contains approx- 
imately a hundred sectors. In principle, this is not a problem, but it can definitely 
become an important issue in practice. The only solution to this issue that we can see 
is a better choice of phase-space parameterization, which leads to a smaller number of 
sectors. 
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